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We  deal  with  the  construction  of  designs  for  linear  multi response 
models.    Three  separate  topics  in  this  area  are  considered. 

The  first  topic  is  the  generation  of  multi response  D-optimal 
designs  when  the  variance-covariance  matrix  of  the  random  error  vector 
is  not  known.    If  E( )  =  F1  ^ ,  1  =  1,  2,  r  are  the  fitted  models 

in  a  multi response  experiment  with  r  responses,  then  a  D-optimal  design 
consists  of  a  set  of  design  points  xj,  _x_2»  •••»  *N  which  maximizes  the 
determinant  of  F' (Z_18IN)F/N  where  F  is  the  block  diagonal  matrix  of 
the  form  F  *  diag(Fls   F2,  Fp)  and  z  is  the  variance-covariance 

matrix  of  the  random  error  vector.  The  design  points  are  chosen 
iteratively  in  a  process  which  involves  the  estimation  of  z  by  a 
consistent  estimator  at  every  iteration. 

The  second  topic  concerns  the  development  of  two  design  criteria  to 
improve  the  power  of  the  multi response  lack  of  fit  test  associated  with 


a  linear  multi response  model.  These  design  criteria  are  multivariate 
extensions  of  the  Aj  and  A2-optimality  criteria  for  the  single  response 
case  discussed  by  E.  R.  Jones  and  T.  J.  Mitchell  (Design  criteria  for 
detecting  model  inadequacy,  Biometrika,  65  [1978],  pp.  541-551).  A 
procedure  is  presented  for  the  sequential  generation  of  a  design  based 
on  the  A2~criterion. 

The  third  topic  is  the  construction  of  robust  designs  to  reduce  the 
effect  of  nonnormality  of  the  error  distribution  on  tests  of  hypotheses 
associated  with  linear  combinations  of  the  parameter  vectors  in  a  linear 
multi response  model. 

Each  of  the  above  topics  is  important  in  the  area  of  multi response 
surface  designs.  Furthermore,  the  techniques  we  develop  could  be  used 
for  better  design  of  multi response  experiments. 


CHAPTER  ONE 
INTRODUCTION 


1.1  A  Note  on  Multi response  Surface  Methodology 
An  experiment  in  two  or  more  responses  are  measured  for  each 
setting  of  a  group  of  independent  variables  is  called  a  multi response 
experiment.  As  an  example,  an  industrial  engineer  may  want  to  study  the 
influence  of  cutting  speed  and  depth  of  cut  on  the  life  of  a  tool  and 
the  rate  at  which  it  loses  metal.  Similarly,  a  medical  researcher 
studying  the  effects  of  complexing  agents  on  the  yield  of  a  certain 
antibiotic  may  also  be  interested  in  the  production  cost,  and  a  food 
scientist  may  wish  to  determine  how  the  flakiness,  gumminess,  and 
specific  volume  of  a  pie  crust  vary  according  to  the  content  of  water, 
flour  and  shortening.  The  settings  of  the  independent  variables  are 
usually  controlled  by  the  experimenter  and  these  variables  are  therefore 
called  controllable  variables. 

Perhaps  the  most  important  objective  of  the  analysis  of  such  data 
is  the  estimation  of  the  relationships  between  the  responses  themselves 
as  well  as  between  the  responses  and  the  controllable  variables  within 
some  specified  region  of  interest.  If  there  are  r  responses  yj,  y2, 
....    yr    and    k    controllable    variables    5lf    K2,  Ck,  these 

relationships,  or  models,  can  be  written  as 
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yui  =  W  V  +  eui'     i=1'  2»         r;     u=1'  2»  N 

(1.1) 

where  yul-  is  the  uth  observation  on  the  ith  response,  jjy  =  (Cul,  Cu2» 
5^)'  is  the  vector  of  controllable  variables  at  the  utn  setting, 
is    the    vector   of   unknown    parameters    associated   with    the  itn 
response,  and  eui  is  the  experimental  error  in  the  utn  observation  on 
the  itn  response  (1=1,  2,  r;  u=l,  2,  N).    In  the  development 

which  follows  we  will  find  it  convenient  to  work  with  unitless  variables 
instead  of  the  independent  (previously  called  controllable)  variables. 
We  define  the  unitless  variables  in  coded  form  as 

Cui"  fi 

xui  =  t —   i=1'  2»  •••»  r;  u=1»  2'  •••»  N 

where  I.  is  the  average  of  Cui  and  Si  is  some  scale  constant. 
Hereafter  we  refer  to  the  xui  as  controllable  variables. 

A  vector-valued  function  consisting  of  the  r  responses,  yj,  y2, 
yr,  is  called  a  multi response  function,  and  the  univariate  models 
given  in  (1.1)  make  up  the  corresponding  multi response  model.  The  exact 
forms  of  the  n/s  are  not  usually  known,  and  will  be  approximated  by 
polynomials  which  are  chosen  by  the  experimenter  in  advance. 
Observations  on  all  r  responses  can  then  be  collected  at  some  specific 
values  of  x_ =  (x1}  x2,  xk)'  determined  by  a  suitable  design.  These 
values  of  x_  are  called  design  points  and  the  choice  of  the  design  is 
based  on  an  appropriate  criterion  which  depends  on  the  model  given  in 
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(1.1).  After  the  data  have  been  collected,  the  unknown  parameter 
vectors  3^-,  i=l,  2,  r  are  estimated  using  a  technique  such  as  least 
squares  or  maximum  likelihood  estimation  and  the  fitted  response  models 
become  completely  determined.    These  models  can  be  used  to 

1.  predict   the   response   values   at   any   given   setting  of  the 
controllable  variables  within  the  experimental  region 

2.  determine  the  optimum  conditions  on  the  set  of  controllable 
variables  that  optimize  the  multi response  function 

3.  test     hypotheses     concerning     the     multi response  model's 
parameters. 

The  predicted  response  functions  can  provide  a  good  representation 
of  the  behavior  of  the  true  responses  depending  on  the  adequacy  of  the 
fitted  model  forms  and  the  precision  with  which  the  parameters  are 
estimated.  Thus  in  most  situations  the  design  criteria  are  intended  to 
improve  the  precision  of  the  parameter  estimates  and  to  increase  the 
power  of  tests  associated  with  detecting  model  inadequacies. 

1.2  Literature  Review 
Aitken  (1935)  suggested  that  the  method  of  generalized  least 
squares  could  be  applied  to  estimate  the  parameters  associated  with  a 
multiresponse  model.  Aitken's  estimators,  however,  depend  on  the 
variance-covariance  matrix  £  associated  with  the  error  vector  which  is 
usually  unknown.  Zellner  (1962)  obtained  estimates  for  the  elements  of 
E  based  on  the  residuals  derived  from  an  equation-by-equation 
application  of  ordinary  least  squares.     This  estimation  procedure  is 


briefly  explained  in  Chapter  Two  and  it  was  shown  by  Zellner  that  the 
parameter  estimates  thus  obtained  were  asymptotically  more  efficient 
than  single  equation  least  squares  estimators.  A  more  general  method  of 
multi response  parameter  estimation  which  applies  to  nonlinear  as  well  as 
linear  models  was  developed  by  Box  and  Draper  (1965),  by  following  a 
Bayesian  approach  using  non-informative  prior  distributions  for  all  the 
parameters  in  the  multiresponse  model,  including  the  elements  of  E.  The 
estimates  for  £  =  (0j,  J^,  &r)' ,  were  obtained  by  minimizing  the 

determinant  of  the  rxr  matrix  V  =  WW  with  respect  to  J3,  where 

W  =  [w  .]  =  [y  .  -  n  .] 
uiJ     L,/ui  uiJ 

and 

^ui  =  VV  -i^'       1  "  1.  2,  ....  r;     u  =  1,  2,  ...  N  . 

The  matrix  V  is  called  the  dispersion  matrix.  This  estimation  procedure 
does  not  require  knowledge  of  £,  thus  eliminating  the  need  for  its 
estimation. 

It  was,  however,  shown  by  Box  et  al .  (1973)  that  this  estimation 
procedure  can  be  adversely  affected  by  the  presence  of  linear 
dependencies  in  the  data.  Therefore  it  is  important  that  such 
dependencies  be  removed  prior  to  the  maximization  of  the  determinant  of 
V.  In  some  cases  the  experimenter  may  not  even  be  aware  of  these 
dependencies.  Box  et  al .  (1973)  gave  a  procedure  to  detect  these 
relationships  among  the  data  by  means  of  an  eigenvalue-eigenvector 
analysis  based  on  the  matrix  DD'  where 


» 

(1.2) 


D"  =  [dui]  =  [yui  -  y.],      i  =  1,  2,  r;     u  =  1,  2,  N 


McLean  et  al .  (1979)  indicated  that  the  presence  of  singularities  in  V 
depended  not  only  on  the  data  but  also  on  the  form  of  the  model.  In 
other  words,  inadequacies  in  the  fitted  models  may  also  jeopardize  the 
technique  of  parameter  estimation.  It  is  thus  quite  desirable  to 
perform  a  test  of  multi response  lack  of  fit  to  detect  any  inadequacies 
in  the  fitted  models.  Khuri  (1983)  developed  a  multi response  lack  of 
fit  test  and  also  a  procedure  that  can  be  used  to  determine  those 
responses  which  are  most  contributors  to  lack  of  fit.  A  brief 
discussion  of  this  test  is  given  in  Chapter  Three. 

In  past  statistical  literature  the  topic  of  multi response  designs 
has  been  largely  ignored  except  for  the  work  done  by  Fedorov  (1972, 
Ch.  5)  on  multi response  D-optimal  designs  when  the  variance-covariance 
matrix  of  the  error  vector  is  known. 

1.3  Objectives  of  the  Present  Work 
The  main  purpose  of  this  research  is  to  construct  designs  for 
linear  multi response  models.  Three  separate  topics  in  this  area  are 
considered.  In  Chapter  Two,  a  sequential  procedure  is  developed  to 
construct  D-optimal  designs  when  the  variance-covariance  matrix  of  the 
error  vector  is  not  known.  In  addition,  a  convergence  proof  for  this 
procedure  is  given.  In  Chapter  Three,  two  design  criteria  are  developed 
to  improve  the  power  of  the  multi response  lack  of  fit  test.  These  are 
the   multivariate    extensions    of   the     A.  and  A?     optimality  criteria 
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defined  by  Jones  and  Mitchell  (1978).  An  iterative  design  algorithm 
which  converges  to  a  A,,-optimal  design  is  also  presented.  In  Chapter 
Four,  the  construction  of  designs  to  reduce  the  effect  of  nonnormality 
of  the  error  distribution  on  tests  of  hypotheses  (robust  multi response 
designs)  is  discussed.  These  tests  deal  with  linear  combinations  of  the 
parameter  vectors  and  the  design  criterion  is  to  minimize  a  multivariate 
measure  of  kurtosis  for  the  controllable  variables,  derived  by  Mardia 
(1971).  The  final  chapter  consists  of  concluding  remarks  and 
suggestions  for  future  research. 


CHAPTER  TWO 
D-OPTIMAL  DESIGNS 

2.1  Introduction 

There  has  been  a  great  deal  of  attention  focussed  on  the  topic  of 
D-optimal  designs,  especially  in  single  response  experiments.  The  main 
importance  of  such  designs  is  that  they  minimize  the  volume  of  the 
confidence  region  of  the  regression  parameter  vector  and  thus  increase 
the  precision  of  the  parameter  estimates.  St.  John  and  Draper  (1975) 
have  written  an  extensive  survey  paper  on  this  subject  that  traces  its 
historic  development  and  contains  a  comprehensive  bibliography.  The 
paper  also  discusses  algorithms  to  construct  D-optimal  designs  including 
the  well  known  ones  developed  by  Fedorov  (1972)  and  Wynn  (1970). 

Construction  of  D-optimal  designs  for  multi response  experiments  was 
considered  by  Fedorov  (1972).  However,  his  algorithm  required  that  the 
variance-covariance  matrix,  £,  of  the  error  vector  be  known.  This  is 
rarely  the  case  in  practice.  The  primary  objective  of  this  chapter  is 
to  develop  an  algorithm  to  construct  multi response  D-optimal  designs 
when  Z  is  not  known. 
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2.2    Model  and  Notation 
In  this  section  we  introduce  our  notation  for  this  chapter.  We 
consider  a  situation  in  which  observations  are  made  on  r  responses  at  N 
experimental    runs    (not   necessarily   all    distinct)    such   that  the  ith 
response  at  the  ut'1  experimental  run  obeys  the  model 

yui  =  VV  -1^  +  eui'    1  "  1»  2'  •••»  r;    u  a  l*  Z>  •••»  N 

(2.1) 

Here  yu^   is  the  observation  on  the  ith  response  at  the  utn  setting, 

Xy  =  (xul»  xu2'  xuk^'  ^s  the  vector  °f  controllable  variables  at 

the  uth  setting  and  belongs  to  the  experimental  region  x;  a  compact 
subset  in  the  k-dimensional  Euclidean  space,      ,  1  =  1,  2,  r  is  a 

vector  of  p^  unknown  parameters  associated  with  the  ith  response  model, 
and  eui  is  the  random  error  in  the  ith  response  for  the  uth  run.  The 
set  of  points  x_lt  x^,  forms  an  N-point  discrete  design  and  is 

denoted  by  DN.  Assuming  that  every  model  in  (2.1)  is  linear  in  the 
parameters,  we  can  write  (2.1)  as 

*ui  =li(xl|)ii  +  eui.    i  =  1,  2,  r;    u  =  1,  2,  N 

(2.2) 

where  fj(jc)  is  a  pi  x  1  vector  which  depends  on  the  form  of  the  response 
function  assumed  by  the  ith  model  and  whose  elements  are  known  functions 
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of  the  controllable  variables  assumed  continuous  within  the  experimental 
region  x.    Using  matrix  notation  the  above  model  can  be  written  as 


Ij  =  Fiii  +  ei      1  -  1,  2,  ....  P 


(2.3) 


where  Yj  is  an  Nxl  vector  of  observations  on  the  itn  response,  F.j  is  an 
Nxp.j  matrix  of  f ul 1  column  rank  with  the  row  containing  the  elements 
of    fj(x  ),    u  =  1,  2,  N,  and       is  an  Nxl  vector  of  random  errors 

associated  with  the  ith  response  (i  =  1,  2,  r) ,  with  E(_e.j)  =  0  and 

E(£jei)  =  a.jj  IN;  i,  j  =  1,  2,  r.    This  implies  that  the  errors 

associated  with  any  single  response  are  homoscedastic  and 
uncorrected.  Thus  the  variance-covariance  matrix  of  the  Nrxl  vector 
£  =  (ej,  £2>  •••»  fr^'  is  9'iven  by 


Var(e)  = 


"l^N  a12TN 
a21!N  °22TN 


arirN  °r2!N 


rr  N 


°11  a12  ' 
a21  °22  ' 


arl  ar2  - 


'lr 


'2r 


rr 


81, 


(2.4) 


=  Z  81 


N 


where  IN  is  the  identity  matrix  of  order  NxN  and  '8'  denotes  the  direct 
product.    The  system  of  equations  in  (2.3)  may  also  be  written  as 
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L 

F.  0  ... 

1 

0 

-1 

-2 

0  F«... 

0 

—2 

-2 

• 
• 

:  ■+ 

• 

• 

• 
• 

0 

F 

*r 

e 

P 

— r 

—  — * 

or 

where   Y'    =  [YJ,  Y£,  r],  B'  =  6^],     e'    =  e£, 

•  ♦.»  £^.]»  and  FD    represents  the  block  diagonal  matrix  on  the  right  hand 
N 

side  of  (2.5).    Let  p  be  the  total  number  of  unknown  parameters  in  the 

r 

system.     Then    p  =  I    p.    is  the  dimension  of  B.    Since  each  F.-  is  of 

i=l  1 

full   column   rank,  FD    is  also  of  full   column   rank  equal   to  p.  As 

N 

suggested  by  Aitken  (1935)  the  technique  of  generalized  least  squares 
can  be  employed  to  obtain  an  estimate  for  _B  given  by 

B  =  [F'  (S^aijF    r!F'  (Z_18I  )Y      ,  (2.6) 
uN  N         N  N  ~ 

and  the  variance-covariance  matrix  of    B  is 
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Var(3)  =  [FJ  (S^H^Fp  ]" 1  . 


(2.7) 


Then  at  point     the  predicted  response  vector  is 


yx(x) 


yP(x) 


o     ...  o 


f£(x)    ...  0 


Ll(x) 


-*rxp 


(2.8) 


or 


y_(x)  =  ♦'(x)e 


with  variance-covariance  matrix 


Varfi(x)]  =  ♦,(x)CFJn(E"1Hn)Fd  r^fx)  (2.9) 


where  y_(x)  =  (y^x),  y2(x),  yp(x))'  and  ♦,(x)  represents  the  rxp 
block  diagonal  matrix  on  the  right  hand  side  of  (2.8). 
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2.3    Design  Theory 
The  concept  of  a  design  measure  defined  below  is  important  to  the 
theory  of  optimal  designs. 

Definition  2.1 

A  design  measure  defined  on  an  experimental  region  x  is  a 
probability  measure,  c(_x)  on  x»  which  satisfies 

C(x)  >  0    and    /  dc(x)  =1     >     x  e  x  (2.10) 
X 

(see  Wynn  (1970),  Kiefer  (1959),  Fedorov  (1972)). 

Given  a  measure  c(x_)  on  x»  the  set  of  points  _x  e  x  at  which 
c(x_)  >  0  is  called  the  support  or  spectrum  of  s(x_).  An  example  of  a 
design  measure  is  as  follows.     Let  _x_i»  2<2,  z  x  and  let  c(x_) 

defined  as 

0        21  *  iLf         1*1.2,...,  s 

C(x)  = 

x^       2L  =  Xj    '     1*1.  2,  ....  s 

s 

where       E    X.  =  l    with    0  <  X.  <  1.    The  design  measure  e(x)  is  said 

1-1    1  1  ~ 

to  be  discrete  if  X^   is  a  rational  number  for  all   i  =  1,  2,  s, 

otherwise,  if  at  least  one  x^  is  irrational,  then  it  is  continuous.  In 

particular,  we  define  the  discrete  design  measure  sD  (x)  by 
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0         **  u  ■  1,  2,  s 

5D  (x)  = 
N  n 

2L=^y    ♦    u  =  l»2,  . . . ,  s 


s 

where  nu  is  the  number  of  replications  at  the  point  Xy  and     z    n    =  N. 

u=l  U 

In  this  manner  we  can  associate  a  design  measure  with  any  given  design. 
For  a  given  design  measure  c(x_)  on  x  and  a  given  variance-covariance 
matrix  Z,  we  define  the  moment  matrix  M(c,  I)  as 


M(C,  I)  =  /  ♦  (x)E_V(x)dc(x)      •  (2.11) 
X 


In  particular,  for  a  discrete  design  measure    cn  (x) 

MUD  .  «)  !    ♦(xu)r"1*'(xu)  (2.12) 

N  u=l 


or  equivalently 


FD  „(r"1,VF0L1 
M(V  E>  ■  —  R  51     •  (2.13) 


2.3.1    Properties  of  M(g,  S) 

For  proofs  see  Fedorov  (1972,  p.  210). 
Let  H  be  the  class  of  all  design  measures  on  x-  Then 
1.     For  ?  s  H,  M(c,  2)  is  positive  semidefinite. 
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2.  The  set  A(E)  of  all  matrices  M(c,  z),  s  e  H  is  convex  i.e., 
if  0  <  o  <  1  and  Hp  e  A  (Z),  then  M  =  (1  -  a)M1  + 
aM2  e  A(Z). 

3.  For  c  e  H,  M(c,  2)  can  be  represented  as 

M(c,  Z)  =    Z    X  *(x  )Z"V(x  )      ,  (2.14) 
u=1    u   -u  -u 

s 

where     1  <  s  <  p\  0  <  X    <  1    with       Z    X    =  1,     and  p'  = 

u  u=l  u 

(p(p  +  l)/2)  +  1. 

In  view  of  (3)  it  is  clear  that  given  any  design  measure 
?,  there  exists  a  design  measure  s'  which  has  a  finite 
support  with  at  most  p'  distinct  points  such  that  M(c,  z)  = 
M(c',  Z).  However  the  fact  that  s'  has  finite  support  does 
not  mean  that  c'  is  associated  with  a  discrete  design,  since 
the  value  of  c'  (x_)  assigned  to  at  least  one  of  these  points 
could  well  be  irrational. 


2.3.2  D-optimality 

Under  the  assumptions  of  normality  of  errors  and  linearity  of  the 
fitted  models,  a  confidence  ellipsoid  for  _S,  of  a  given  confidence 
coefficient  has  the  form 


{0:  (g_  -  6)'M(?,  Z)(j3  -  j3)  <  constant} 
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where  e_  is  the  generalized  least  squares  estimate  of  _B  (see  Silvey 
[1980,  p.  10]  for  the  single  response  case).  A  more  precise  estimate 
for  S_  is  obtained  by  making  this  ellipsoid  as  small  as  possible.  Since 
the  volume  of  the  ellipsoid  is  proportional  to  |M(c,  E)|-1/2,  this  can 
be  achieved  by  choosing  a  design  measure  which  maximizes  |M(c,  E)|. 
Such  a  design  measure  is  said  to  be  D-optimal  and  can  be  defined  as 
follows. 

Definition  2.2 

Let  H  be  the  class  of  all  design  measures  on  x«  Then  c*  is  called 
D-optimal  with  respect  to  E  if 

|M(c*.  Z)|  =  sup|M(C,  I)|      .  (2.15) 


2.4  Construction  of  D-optimal  Designs  When  E  is  Known 
Fedorov's  algorithm  to  construct  D-optimal  designs  (when  E  is 
known)  is  a  sequential  procedure,  in  the  sense  that  design  points  are 
chosen  one  at  a  time.  It  is  based  on  a  certain  equivalence  theorem.  In 
this  section  we  include  the  sequential  procedure  along  with  the 
statements  of  the  equivalence  theorem  (Theorem  2.1)  and  the  theorem 
concerning  convergence  of  the  sequential  procedure.  The  interested 
reader  is  referred  to  Fedorov  (1972,  Ch.  5)  for  a  more  detailed 
discussion. 
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Theorem  2.1    (Equivalence  Theorem)  (Fedorov,  1972,  p.  212) 
Let  H  be  the  class  of  all  design  measures  on  x»  and  let 

V(x,  C,  E)  =  ♦,(x)M"1U,  r)«(.(x)      ,      c  e  H     .  (2.16) 

The  following  conditions  are  equivalent. 

1.  The  design  measure  c*  is  D-optimal . 

2.  sup  trCE'Vx,  C*,  Z)]  =  inf  sup  trCl'Mx,  ?,  E)]  . 
x.ex  CeH  xex 

3.  sup  trCi^v^,  <;*,  £)]  =  p. 

X^X 

2.4.1    The  Sequential  Procedure  When  E  is  Known  (Fedorov,  1972,  Ch.  5) 

For  convenience,  throughout  the  rest  of  this  chapter  we  will  write 
5N  ^or  5DN*  ^s  described  by  Fedorov  we  start  with  a  nondegenerate 
design  DNq,  i.e.,  m(?Nq»  z)  is  nonsingular.  By  successive  addition  of 
points  to  DNq  we  will  construct  a  sequence  of  designs  dNq+1»  dNq+2»  ••• 
which  converges  to  a  D-optimal  design.  The  steps  in  this  procedure  are 
as  follows. 

1.  Find  x         in  x  such  that 

trCE"1V(xfJ        c    ,  1)1  =  sup  tr[E-1V(x,  r    ,  s)]  . 
0  0  xex  0 

2.  Form  the  new  design  DN  +j  by  adding       +j  to  DN  . 
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3.  Continue  this  procedure  and  obtain  the  nested  sequence  of 
designs  D^  DNQ+1  "*  wnere  DN  is  obtained  from  DN-1  by 
adding  the  point  _x^  in  x  which  satisfies 


trCl'Vx^,  C„  v  1)1  =  sup  tr[E-1V(x,  s        E)]      .  (2.17) 

The  stopping  point  for  this  procedure  is  reached  when  we  find  x^,  in  x 
which  satisfies 


trCE^V^.  ,  zn,_v  E)]  -  p  <  6 


where  6  is  a  small  positive  number  chosen  a  priori. 

Theorem  2.2  (Fedorov,  1972,  p. 216) 

The  iterative  procedure  described  above  converges  to  a  D-optimal 
design,  i.e., 


lim|M(5N,  E)|  =  |M(?*,  E)| 


where  s*  is  a  D-optimal  design  measure. 


2.5    Construction  of  D-optimal  Designs  When  E  is  Not  Known 

As  was  clearly  seen  the  moment  matrix  M(c,  E)  does  depend  on  E. 

For   this   reason   the  above  procedure  cannot  be  applied  unless   E  is 

known.     If  E  is  unknown,  a  consistent  estimator  of  E  can  be  used  to 
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construct  a  sequence  of  random  design  measures  which  converges  in 
probability  to  a  D-optimal  design  measure  with  respect  to  E.  (For 
definition  of  convergence  in  probability  see  Appendix  A).  Such  a 
consistent  estimator  (see  Appendix  C  for  proof)  was  proposed  by  Zellner 
(1962)  and  is  given  by    E^  =  [o?.]  where 

*  U,  -        ) ' (Xj  "  Fjlj)      ,     U  =  1.  2,  ....  r, 

(2.18) 

~  -1 

and  =  (Fi!Fi )  FjY^  is  the  usual  single  equation  least  squares 
estimator,  i  =  1,  2,  r.    Now  suppose 

A  =  (diag  E-V1/2  E-^diag  e"1)"172  (2.19) 

where  diag(P)  denotes  the  diagonal  matrix,  whose  diagonal  elements  are 
equal  to  the  diagonal  elements  of  the  square  matrix  P.  Then  the 
following  theorem  proves  that  a  0-optimal  design  with  respect  to  E  is 
equivalent  to  a  D-optimal  design  with  respect  to  A. 

Theorem  2.3 

|M(c,  E)|  =    n    (<T11)Pl|M(c,  A"1)!,  C  e  H 
i  =1 

where  E"1  =  [aij]  and  p.  is  the  number  of  parameters  in  the  model  Y-  = 

f  — ] 

^■fjLj  +  L\ •     Tne  proof  is  given  in  Appendix  D. 
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It  can  be  seen  from  the  expression  on  the  right  hand  side  that  only 
lMU»  A"1) |  depends  on  the  design  measure  c  Hence  a  design  measure 
which  maximizes  |M(c,  A"1)!  will  essentially  maximize  |M(c,  Z)|.  It  can 
be  also  shown  that  the  results  stated  in  Section  2.4  hold  true  when  Z"1 
is  replaced  by  A  and  thus  an  asymptotically  D-optimal  design  measure 
with  respect  to  A  can  be  obtained  using  the  sequential  procedure 
described  in  Subsection  2.4.1.  However,  when  Z  is  not  known,  A  is  also 
not  known  and  so  we  modify  this  procedure  so  that  A  is  replaced  by 

AN  =  (diag  1/2  i"1  (diag  i^)"1/2  (2.20) 

which  will  be  updated  every  time  a  new  design  point  is  chosen.  This 
procedure  is  described  in  Section  2.6. 

The  construction  of  a  D-optimal  design  with  respect  to  A  (instead 
of         is  more  desirable  for  two  reasons. 

1.  Since  the  elements  of  A  lie  between  -1  and  1,  A^  reaches 
stability  much  faster  than  z^  thereby  giving  rise  to  a 
rapidly  converging  point  in  the  sequential  procedure 
described  below. 

2.  Due  to  the  fact  that  the  diagonal  elements  of  A  are  all  equal 
to  1  the  number  of  elements  of  A  to  be  estimated  is  reduced 
by  r. 

2.5.1    The  Sequential  Procedure  When  Z  is  Not  Known 

1.     Start  with  an  initial  nondegenerate  design  D^   with  Nq  points 

and  compute    z       and    Aw      (based  on  observations  on  all  r 
n0  N0 
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responses  measured  at  all  the  points  of  D^)  using  (2.18)  and 
(2.20). 

2.     Construct    DN    .     by  adding    Xj.        to    Dw     where  xw 
satisfies 


tr[AN  V(x^         c    ,  A"1)-  «  sup  tr[A    V(x,  ?    ,  A"1)-  , 
"0    ^0  A      M0      N0        xex        N0  N0  N0 


and  obtain    cN  +1    by  assigning  mass  1/(NQ+1)  to  each  design 
a  0 

point  in     DK.  ±. 

NQ+1 . 

3.  Once    D^+^,  N  >  1,    is  obtained,  compute    ZN  +N    and    A^  +N 
(based  on  the  observations  on  all  r  responses  measured  at  all 

the  points  of    DN    N)    using  (2.18)  and  (2.20). 
.  0 

4.  Construct     D^+n+l     by    adding     XflQ+N+1     to     DN    N  which 
satisfies 


tr^AN0+NV(^N0+N+l'  ?NQ+N'  ^q+N5' 


=  sup  tr^VU,  1%+H,  A-J+N). 


and  obtain    5N  +N+1,    by  assigning  mass  1/(NQ+N+1)  to  all  the 

0  a 
design  points  in 

The  stopping  rule  for  this  procedure  is 

sup  tr[AN, V(x,  c„, ,  A"})-  -  p  <  6 
x_ex 

o  is  chosen  beforehand  and  N '  is  some  positive  integer  >  Nn. 
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We  now  consider  convergence  of  the  procedure  given  above.    Let  A  be 

the  set  of  all  rxr  symmetric  matrices  A  =  Ca-jj]  such  that  a^  ■  1,  1  ■ 

1,  2,  r  and    -1  <  ai  .  <  1,        1  <  i  <  j  <  r.    Define  _b(A)  =  (bj, 

b2,  bpi)'  where  r'  =  r(r-l)/2,  to  be  the  vector  which  consists  of 

the  elements  of  A  above  its  diagonal  taken  in  order  from  left  to  right 

for  each  row  starting  with  the  first.    We  call  b_(A)  the  r' -dimensional 

vector  associated  with  A.    Note  that  since  AeA  is  symmetric  and  all  its 

diagonal  elements  are  equal  to  1,  A  can  be  completely  described  by  the 

r'  elements  a^,    1  <  i  <  j  <  r.     It  is  clear  that  A  and    {A^,  N  >  1} 

defined  in  (2.19)  and  (2.20),  respectively,  belong  to  A.     Also  it  is 

shown    in    Appendix    E    that     AN     is    a    consistent    estimator    of  A. 

Equivalents,   if    eN    is   the  Euclidean  distance  between    b(AJ  and 

b_(A),    then    eN    converges  in  probability  to  0. 

It  is  conjectured  that  the  convergence  in  probability  to  0  of  e^ 

implies  that  the  above  sequential  procedure  converges.    However,  we  are 

able  to  prove  convergence  of  the  sequential  procedure  only  by  assuming  a 

N  « 

stronger  condition  that     Z    e     converges  in  probability  to  some  random 

u  =  l  u 

variable  e.    This  can  be  formally  stated  as  follows. 
Theorem  2.4 

N  . 

Suppose    that       I    e      converges    in   probability   to   some  random 
u=l 

variable  e.  Then  for  a  given  6  >  0,  there  exists  an  integer  N  >  0  such 
that 

sup  tr[A  V(x,  CN,  A"J)«  -  p  <  6 

xex 
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with  probability  1  (see  Appendix  A  for  definition  of  convergence  with 
probability  1). 

The  proof  is  given  in  Appendix  F. 

2.6    A  Numerical  Example 
We  now  present  an  example  to  illustrate  the  procedure  described  in 
Section  2.6.    We  consider  an  experiment  with  two  responses  which  depend 
on  three  controllable  variables.    The  experimental  region  x  =  (x  3  Uj, 
x2,  x3):  |x.j|  <  1.73.  i  =1,  2,  3}  and  the  fitted  models  are 

yl  =  B10  +  311X1  +  e!2X2  +  e!3x3  +  *112xlx2  +  eH3X!X3  + 


2  2 
0111X1  +  3133x3  +  el 


(2.21) 


y2  =  S20  +  S21xl  +  822x2  +  3212xlx2  +  B211xl  +  B222x2  +  £2 


Thus  £j  and  ^  are  of  dimensions  pj  =  8  and  p2  =  6  and  the  total  number 
of  unknown  parameters  is  p  =  14.  The  matrices  Fj,  F2  consist  of  eight 
and  six  columns,  respectively,  and  the  matrix  Fn    is  of  the  form 


1 


0 


i he  error  vector  e  =  (ej ,  e£ ) '  is  assumed  to  be  normally  distributed 
with  zero  mean  vector  and  variance-covariance  matrix    Z8I  where 
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E  = 


2  .4 
.4  1 


The  matrix  I  will  be  used  to  generate  the  normal  error  values  and  e2 
with  the  stated  variance-covariance  matrix  but  will  not  be  used  in  the 
sequential  generation  of  the  D-optimal  design.  Also 


.-1 


,543 


.217 


.217 


1.086 


and 


A  = 


1 

•.28 


■.28 
1 


The  error  vector  £  was  generated  by  means  of  a  SAS  procedure  which 
produces  normally  distributed  random  vectors.  The  initial  design 
consisted  of  NQ  =  15  points  and  is  given  in  Table  2.1.     At  the  Nth 


iteration  a.: 


»  1  <  1,  j  <  2   was  calculated  by 


(2.22) 


where  Ri  =  F^FjF.)"1  F! ,    i  =  1,  2,  and  _Ek,  Fk,  k  -  1,  2,  consists  of 

(Nq+N)   rows.      Observe  that   in  situations  where   response  values  are 
,N0+N 

available,  a..        should   be   calculated   using    (2.18).      Using  (2.20) 
a12     was  then  calculated  and  these  values  are  given  in  Table  2.2.  The 
maximization    at  each  iteration  was  carried  out  using  a  computer  program 
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Table  2.1 

The  Initial  Design 

Xj                           x2  x3 

1                   1  -1 

1                   1  1 

1                  -1  -1 

1                  -1  1 

-1                  1  -1 

-1                  1  1 

-1                  -1  -1 

-1                  -1  1 


1.68  0  0 

0  1.68  0 

0  0  0 

0  0  0 

0  0  0 

0  0  0 

0  0  0 
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Table  2.2 

The  Estimates  for  a^  ancl  »  i ,  j = 1 ,  2} 
~NQ+N-1 

-i°2  ^+N-1  ^+N-1  -J02+N-l 


15 

.302 

.752 

.958 

.256 

16 

.334 

.722 

1.03 

.280 

17 

.220 

.680 

1.06 

.186 

18 

.214 

.648 

1.01 

.172 

19 

.024 

.653 

1.03 

.018 

20 

.278 

.657 

1.36 

.266 

21 

.214 

.660 

1.35 

.202 

22 

.256 

.670 

1.30 

.238 

23 

.218 

.660 

1.33 

.204 

24 

.080 

.703 

1.32 

.076 

25 

.090 

.678 

1.29 

.084 

26 

.114 

.652 

1.30 

.104 

27 

.212 

.660 

1.28 

.194 

28 

.298 

.638 

1.29 

.270 

29 

.292 

.616 

1.28 

.260 

30 

.284 

.596 

1.33 

.252 

31 

.298 

.584 

1.40 

.270 

32 

.304 

.623 

1.55 

.300 

33 

.284 

.699 

1.54 

.294 

34 

.332 

.619 

1.89 

.360 

35 

.318 

.632 

1.84 

.344 

36 

.320 

.621 

1.85 

.342 

37 

.332 

.605 

1.83 

.350 
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Table  2.3 


The  Augmented  Design  Points  for 

a  u-uptimai 

uesi gn 

NQ+N 

tr\.+N-lV(-,?Nn+N-l  'ANn+N-l ) 
u              u  u 

xl 

X3 

16 

105.00 

1.724 

1.724 

-1.728 

17 

101.47 

1.729 

1.727 

-1.703 

18 

88.74 

1.728 

-1.729 

-1.720 

19 

84.29 

-1.728 

-1.725 

-1.680 

20 

50.73 

1.729 

1.729 

1.729 

21 

52.00 

-1.725 

-1.723 

1.715 

22 

42.42 

-1.730 

1.721 

1.729 

23 

48.80 

1.730 

-1.729 

1.729 

24 

22.95 

-1.730 

1.730 

.026 

25 

23.95 

1.708 

1.660 

1.692 

26 

23.35 

1.730 

-1.730 

-.045 

27 

25.19 

-1.729 

-1.730 

-1.728 

28 

20.55 

-1.730 

-.096 

1.730 

29 

23.07 

1.729 

1.724 

-1.729 

30 

21.34 

1.730 

-1.730 

.153 

31 

20.58 

-.154 

1.730 

-1.730 

32 

21.14 

1.690 

-1.702 

1.658 

33 

20.48 

1.727 

-1.710 

-1.726 

34 

20.75 

-.101 

-1.730 

1.730 

35 

19.92 

-1.723 

1.601 

1.695 

36 

20.00 

1.729 

1.729 

1.722 

37 

16.56 

1.729 

-1.703 

1.698 
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called  SEARCH  written  by  Conlon  (1979).  The  program  is  based  on  the 
controlled  random  search  procedure  introduced  by  Price  (1977).  This 
procedure  uses  a  random  search  to  locate  an  optimal  point  from  among  a 
collection  of  points,  the  number  of  which  is  predetermined  by  the 
user.  The  augmented  design  points  are  given  in  Table  2.3.  Observe  that 
the  procedure  has  succeeded  in  achieving  a  substantial  reduction  in  the 
' tr '  value,  (105  to  16.56).  Also  the  difference  between  the  final  value 
of  16.56  and  the  anticipated  value  of  14  is  relatively  small.  In 
addition  it  can  be  seen  that  the  values  of  a,2  have  stabilized.  This 
demonstrates  the  effectiveness  of  our  procedure  in  constructing  D- 
optimal  designs  when  £  is  not  known. 

In  Table  2.3  the  last  entry  indicates  a  'tr'  value  of  16.56  at  the 
22nd  iteration.  When  further  iterations  were  carried  out  no  further 
drop  was  noted  in  the  'tr'  value.  This  may  be  due  to  the  fact  that  the 
number  of  points  used  for  the  random  search  may  have  been  insufficient 
to  locate  the  true  optimal  point. 


CHAPTER  THREE 

DESIGN  CRITERIA  TO  INCREASE  THE  POWER  OF  LACK  OF  FIT  TEST 


3.1  'Introduction 
In  the  past,  considerable  attention  has  been  focussed  on  the  topic 
of  model  inadequacy  in  linear  regression.  Its  importance  is  described 
in  Box  and  Draper  (1959).  Once  a  single  response  experiment  is 
performed  assuming  a  certain  model  is  true,  a  test  to  detect  lack  of  fit 
of  the  fitted  model  can  be  contemplated  (see  Draper  and  Smith,  1981, 
Ch.  2).  Authors  such  as  Aitkinson  (1972),  Aitkinson  and  Federov  (1975) 
and  Jones  and  Mitchell  (1978)  have  come  up  with  design  procedures  to 
improve  the  detection  of  model  inadequacy  in  single  response 
experiments. 

Detection  of  model  inadequacy  is  an  important  issue  in  the 
multi response  case  as  it  is  in  the  single  response  case.  Khuri  (1983) 
developed  a  test  for  lack  of  fit  for  a  linear  multiresponse  model.  As 
mentioned  by  Khuri  (1983),  the  inadequacy  of  the  fitted  model  can 
seriously  affect  the  estimation  procedure  of  regression  coefficients 
developed  by  3ox  and  Draper  (1965).  A  detailed  discussion  of  this 
procedure  and  the  problems  associated  with  it  is  given  in  Box  and  Draper 
(1965),    Box  et  al .    (1973)  and  McLean  et  ai  .  (1979).    Khuri  (1983)  also 
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states  that  the  multi response  lack  of  fit  test  provides  a  comprehensive 
assessment  of  the  adequacy  of  the  multi response  model.  It  can  detect 
lack  of  fit  of  the  overall  model  even  in  situations  where  the  usual 
single  response  lack  of  fit  test  performed  on  each  individual  response 
separately  fails  to  reveal  any  inadequacies  in  the  corresponding 
model.  This  fact  was  well  illustrated  by  Khuri,  by  performing  the 
multiresponse  lack  of  fit  test  on  data  obtained  from  a  study  made  by 
Kizer  et  al .  (1978)  on  the  behavior  of  a  fluidized  bed  reactor  for  the 
catalytic  oxidation  of  benzene  to  mallic  anhydride.  Khuri  (1983)  also 
gave  a  procedure  to  determine  which  responses  are  responsible  for  lack 
of  fit  when  the  multiresponse  test  for  lack  of  fit  is  significant. 

The  purpose  of  this  chapter  is  to  develop  design  criteria  to 
increase  the  power  of  the  multiresponse  test  of  lack  of  fit  mentioned 
above.  The  model  assumptions  and  the  notation  in  this  chapter  follow 
Khuri  (1983)  and  are  given  below  for  completeness. 

3.2    Model  and  Notation 
Let  N  be  the  total  number  of  experimental  runs  and  r  be  the  number 

of  responses.    We  assume  that  each  response  is  dependent  on  all  or  some 

of    the    k    controllable    variables    denoted    xls    x2,            xk.  The 

dependencies    of    the    responses    on    the    controllable    variables  are 

expressible  in  terms  of  polynomials  of  possibly  different  degrees.  The 
fitted  io11  response  model  can  be  represented  as 


30 


where  is  an  Nxl  vector  of  observations  on  the  itn  response,  Ea(yj) 
denotes  the  expected  value  of  Y_|  under  the  fitted  model  (the  subscript 
'a'  denotes  assumed),  Xi  is  an  Nxp^  matrix  of  rank  p^  of  known  functions 
of  the  settings  of  the  controllable  variables,  and       is  a  p^xl  vector 

of  unknown  constant  parameters  (i  =  1,  2  r). 

The  model    for  the  true   ith  response  mean   (i  =  1,  2,  r)  is 

assumed  to  be  of  the  same  form  as  the  fitted  model  but  possibly  contains 
terms  in  addition  to  those  in  the  fitted  model.  It  can  be  written  as 
follows. 

W  =  Xi*i  +  Zili,    1  -  1.  2.  ....  r  (3.2) 

where  Et(_YJ-)  denotes  the  expected  value  of  Yj  under  the  true  model,  Zi 
is  an  Nxqi  matrix  of  known  functions  of  the  settings  of  the  controllable 
variables,  and  jj  is  a  vector  of  unknown  constant  parameters.  If  the 
fitted  model  (3.1)  is  correct,  then  r,  will  be  equal  to  the  zero  vector. 

The  equations  given  in  (3.1)  and  (3.2)  can  be  written  to  represent 
the  overall  multi response  model  in  the  following  manner. 

Ea<Y)  -  XB  (3.3) 


Et(Y)  =  XB  +  zr 


(3.4) 
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where  Y  =  [Y^  ^  ...:  Y^,  X  =  [X^  ^:  Xp],  Z  =  [Z^  1^. 

Zp],  B  =    DiagD^,  B2,  ....  j^J,    and    r  =  DiagC^,  t2.  •      Yp].  The 

matrices    Y,    X,    Z,    and    r   are   of   orders    Nxr,    Nxp,    Nxq,    and  qxr, 

r  r 
respectively,   where     p  =    z  p.  <  N,      q  =    z    q  ,     and   X  is   of  rank 

i=l  1  i=l  1 

p(<  p).    The  rows  of  Y  are  independent  observations  from  multivariate 

normal  populations  with  a  common  nonsingular  variance-covariance  matrix 
I  of  order  rxr.  Under  the  true  model,  Y  has  a  mean  given  by  (3.4),  and 
is  referred  to  as  the  normal  data  matrix  with  mean  XB  +  Zr  and  variance- 
covariance  matrix  I^BI.    This  is  written  symbolically  as 

Y  ~  N(XB  +  zr,  I  II) 

In  the  development  of  the  lack  of  fit  test  it  is  assumed  that 

replicated  observations  are  available  on  all  r  responses.    Without  loss 

of  generality  it  will  be  assumed  that  such  replicated  observations  are 

obtained  at  each  of  the  first  n  design  points,  where    1  <  n  <  N.  The 

number  of  repeated  observations  at  the  ith  design  point  is  denoted  by  v- 

n 

(v.  >  2,  i  =  1,  2,  ....  n).     Let    v  =    Z    v.  (v  <  N).    The  matrices  Y, 

i-1  1 

X,  and  Z  in  (3.3)  and  (3.4)  are  then  partitioned  as 


"Y(D" 

>>' 

*Z(D" 

Y  = 

,   x  = 

,   z  = 

y(2) 

x(2) 

> 

z(2) 

— 
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where  Yd),  X(D,  and  Z(D  are  matrices  of  orders  vxr,  vxp,  and  vxq, 
respectively,  and  are  associated  with  the  replicated  portion  of  the 
design,  whereas  the  matrices  Y(2),  X(2),  and  Z<2)  are  associated  with 
the  nonrepl icated  portion. 

3.3  Development  of  Khuri's  Lack  of  Fit  Test 
Khuri  (1983)  developed  a  multivariate  lack  of  fit  test  by  first 
transforming  the  normal  data  matrix  Y  into  a  single  response  normal  data 
vector  and  applying  the  usual  single  response  test  for  lack  of  fit 
described  in  Draper  and  Smith  (1981,  Ch.  2).  The  null  hypothesis  tested 
by  the  multivariate  lack  of  fit  test  is 

H0:  [IN  "  VX0X0)_lx0]Zr  =  0  (3-5) 

where  XQ  is  a  Nxp  matrix  of  rank  p  whose  columns  form  a  basis  for  the 
columns  of  X,  which  is  of  rank  p. 

The  null  hypothesis  HQ  is  the  hypothesis  of  zero  lack  of  fit  (or 
adequacy  of  fit)  of  the  model  (3.3)  when  in  reality  the  true  model  is 
given  by  (3.4).  The  test  statistics  to  test  this  null  hypothesis  as 
developed  by  Khuri  (1983)  is  explained  below.  Let 

K  =  Diag(Kr  K£,  Kn ,  0) 

where  K  is  a  block  diagonal  matrix  of  order  NxN,  0  is  a  zero  matrix  of 
order    (N-v)  x  (N-v),  and 
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K.  =  Iv    -  (l/v.)Jv     ,    1=1,  2,  n      .  (3.7) 

i  1 


Here  Iv.  is  the  identity  matrix  of  order    v.  xv.    and  J     is  the  matrix 
of  ones  of  order    v.xv.  (i  =  1,  2,  . ..,  n). 
Define  the  matrices      and  Go  as 


Gl  =  Y'^N  "  "  K]Y  (3.8) 


and 


G2  =  Y'KY     *  (3.9) 


An  appropriate  test  statistic  to  test  HQ  given  in   (3.5)  is  the 
maximum  eigenvalue  of    G.GI1,    denoted  by    em3  (G.G:1),    and  the  test 

A    £  maX       1  c. 

procedure  is  reject  HQ  at  the    a-level  of  significance  if 

W^1)  >  \  (3.10) 

where     Xq       is   the  upper  100a  %  point  of  the  null   distribution  of 

emax^GlG2^  (i,e"  the  distribution  of  e^GjG'1)  under  H0^'  This 
multivariate  test  is  known  as  Roy's  largest  root  test.  Three  other 
multivariate  test  statistics  are  also  available  to  test  HQ.  They  are 
(1)  Wilks'   likelihood  ratio,  W=  |G2[/|G1  +  G2j,    (2)   Pillai's  trace, 
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tr[G1(G1+G2)-i],  and  (3)  Hotel  1  i ng-Lawl ey 's  trace,  tr^G" L)  where 
'I  |'  and  'tr'  respectively  denote  the  determinant  and  the  trace  of  a 
matrix. 

3.4    Development  of  Design  Criteria 

It    is    shown    by    Khuri    (1983)   that   G2  has   the  central  Wishart 

n 

distribution  with    vpfr  =    E    (v.  -  1)    degrees  of  freedom;  also  that  Gi 

i=l     1  1 
is  independent  of  G2  and  has  the  noncentral  Wishart  distribution  with 

vLp  =  (N-p-vpE)  degrees  of  freedom  and  a  noncentral ity  parameter  matrix 
given  by 

°-  r1r'Z'[lN  -  Wo)-\nt    .  o.n) 

It  is  known  that  the  power  of  the  lack  of  fit  test  based  upon  any  of 
these  four  test  statistics  mentioned  above  is  a  montone  increasing 
function  of  the  eigenvalues  of  ft  (see  Roy  et  al.,  1971,  p.  68). 
Therefore,  the  power  of  this  test  can  be  increased  by  increasing  the 
eigenvalues  of  ft.  This  can  in  turn  be  achieved  by  maximizing  the  trace 
of  ft.  This  is  true  since  increasing  the  trace  of  ft  will  increase  at 
least  one  of  its  eigenvalues.  However,  the  choice  of  the  design  which 
maximizes  the  trace  of  ft  depends  on  the  matrices  £  and  r  which  are 
unknown.  Thus  we  are  faced  with  the  problem  of  finding  an  expression 
independent  of  Z  and  r.  whose  maximization  will  result  in  an  increase  in 
the  trace  of  ft.    This  expression  is  found  as  follows. 


35 


From  Theorem  B.3  in  Appendix  B  we  know  that 

tr(n)  =  tr{z-1/2r1z,[iN-xn(x'xn)-1x']zrz"1/2} 


>  wz_1/ >tr{r'z'ivvw~Vr}  •  (3-12) 

Recall  that  Z  =  [Z^  Z2:  . ..:  Zp]  where  the  columns  of  Z^ ,  correspond  to 

the  extra  terms  in  the  true  model  for  the  ith  response  (i  =  1,  2,  .... 

r).    It  is  easy  to  see  that  more  than  one  column  in  Z  can  be  assigned  to 

the  same  term  in  the  extra  portion  of  the  multi response  model  i.e., 

more  than  one  model   can  have  the  same  term.     Let  ZQ  be  the  matrix 

obtained  from  Z  by  deleting  those  columns  which  occur  more  than  once  so 

that  each  term  in  the  true  model  corresponds  to  exactly  one  column  in 

Zq.    For  example  let  us  consider  a  situation  with  two  responses  and  let 

the    columns    of    Zj    and    Z2    correspond    to    terms      {x2,  x2}  and 
2  2 

{xr  x2,  x1x2>,     respectively.     Then  Z  consists  of  five  columns  with 

columns   one  and   three   corresponding   to     x2,     columns   two   and  four 

2 

corresponding  to  x2  and  column  five  corresponding  to  Xjx2  so  that  Zq 
is  obtained  by  deleting  columns  three  and  four. 

In  general,  Z1 ,  i  =  1,  2,  r,  can  be  written  as 


Zi  ~  Z0Hi      '      i  =  1,  2,  r 


(3.13) 
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Here  ZQ  is  of  order  Nxp1  where  pj  is  the  number  of  columns  in  ZQ.  The 
Hi  matrices  are  of  order  p^q^  and  their  elements  consist  of  zeros  and 
ones .  Now 

tr<r'z'tiN  -  Wo)-\nvy 


"HiZi 

II 

=  tr 

H2Z0 

"N-MWlxo]CzoHr-:ZoHr^ 

12 

• 

• 

• 

• 
• 
• 

• 

• 

• 

V 

H'Z' 
r  0 

=  tr 


Y2 


Y*  H 


=  fJ1  Ii,H;AHiIi    where    A  =  ZJd^CX^X^^XJjZ 


CtJHJ:  y£  H£:  .. 


•i;  H;](ir«A) 


Hlll 

H2l2 


HpY 
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where  the  vector    [y'  HI:  yi  HI:  ...:  y'H1]    contains  the  elements  of 

—i    i     c    c  — r  r 

Y-j  >  i  =  1,  2,  r  and  possibly  some  zeros.    We  want  to  move  all  the 

zero  elements  of  [yj  HJ:  j£H£:  ]J  H']  to  the  right  and  the 
nonzero  ones  to  the  left.    This  can  be  accomplished  by  writing 


(>iHi:  l2H2:  iHr]  =  Cy':  °,]C 


(3.14) 


where    y'  =  [y{:  yL:  ...:  yl]       ,    0  is  a   (rprq)xl  vector  of  zeros, 

lxq 

and  C  is  an  orthogonal  matrix  which  has  the  function  of  reordering  the 
elements   of     [y{  HJ:  y£H£:  ...  Y^].     It   is   obtained   by  a  suitable 
rearrangement  of  the  columns  of  the  identity  matrix. 
Therefore 


tr{r«z-[iN  -  x^xjx^xyzr) 


=  [y1 :  0']  C(Ir8A)C 


(rpj-q)xq 


y'L(I  8A)L'y  where  L  =  [I  :  0  ,  ,]C 
-      r        -  q  qxtrp^q)-1 


(3.15) 
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Thus,  inequality  (3.12)  can  be  written  as 


-1/2 


h'  L 


(IpBA)  L'y 


Now  since  emin(Z"  '  )  is  a  constant,  the  above  inequality  implies 
that  maximization  of  the  quantity 


will  result  in  increasing  the  trace  of  Q.  However,  the  choice  of  design 
to  maximize  A  depends  on  y,  which  is  unknown.  In  order  to  overcome  this 
problem  we  apply  the  maximin  method  which  was  proposed  by  Aitkinson  and 
Fedorov  (1975)  and  used  by  Jones  and  Mitchell  (1978)  in  the  single 
response  case.  The  maximin  method  consists  of  choosing  a  design  which 
maximizes  Aj,  the  minimum  of  A  over  a  specified  region  A  in  the 
2-space.  The  specification  of  this  region  A  depends  on  a  quantity  t 
which  can  be  considered  as  a  measure  of  the  inadequacy  of  the  fitted 
model.  Suppose  the  uth  rows  of  Xi  and  Zi ,  (i  =  1,  2,  r; 
u  =  1,  2,  N)  in  (3.2)  can  be  represented  by    f.' (x  )    and    q!(x  ) 

—1  — U  -2-1  v— U  '  ' 

respectively.     Then  the  fitted  and  true  response  functions  associated 


A  =  Y'L(IrBA)L'Y 


with  (3.1)  and  (3.2)  are  f^xjB- 
....  r  respectively.    We  express  t  as 


and    f_\(x)pA  +cLj(xhi, 


i=l,  2, 


t  =  y'  T  y 


where  y 


Cy^  •  Y^2  •    •  •  •  • 


(3.16) 
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and  _ 
T  = 


1 

T, 


T 

r 


,  with  T.=  v\2  -  v'n  (v^rV^,  1  =  1,  2,  ....  r. 


j  -> 

Here  ,  k,l  =  1,  2,  are  the  region  moment  matrices  defined  by  = 
s  /  £i  (2L)£j  (x  )dx ,  mJx  =  S  /  £i  (x)f j  (x)dx,      ujj  =      S  /  f .  (x)f (x)dx, 


x 

and 


u\2  =  S  /  f-(xjg!(x)dx     where    S"1    =    J  dx     and    x  denotes  the 

experimental  region.  This  is  the  multi response  extension  of  the 
expression  for  t  given  by  Jones  and  Mitchell  (1978).  Now  since  T,, 
i  =  1,  2,  r  is  positive  definite,  T  is  also  positive  definite. 

Hence  it  is  clear  that  t  will  be  equal  to  zero  if  and  only  if  y  is  equal 
to  zero.  This  happens  when  the  fitted  model  is  true.  So  t  is  a  measure 
of  the  inadequacy  of  the  fitted  multi response  model  given  in  (3.3)  when 
in  reality  the  true  model  is  (3.4).  It  is  positive  whenever  the  fitted 
model  is  inadequate  and  is  zero  otherwise. 

3.4.1    A i -optimal ity 

Since  t  will  be  positive  if  the  fitted  model  is  inadequate,  then  in 
such  a  case  x  >  6  for  some  constant  5  >  0.  We  define  a  =  {y:  y'Ty  >  <5}. 
The  first  design  criterion  is  to  maximize  Aj  where 


A,  =  inf  y'L(I  »A)L*t  . 
YeA  r  - 


Recall  that    A  =  ZTl    -  X  fx'x 
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This  is  a  multi response  extension  of  the  A^  optimal ity  proposed  by 

Jones  and  Mitchell   (1978)  and  the  minimum  value  of  the  quadratic  form 

y_'l(l  8A)L't     over     A    occurs   on  the  boundary   of  A    which   is  the 

9A 

contour  t  =  6,  since        =  0  only  at  the  origin  which  is  not  an  interior 

point  of  A.  Furthermore,  as  in  Jones  and  Mitchell  (1978),  it  can  be 
shown  that 

Al  =  ^min^V^1']      '  (3.17) 

Since  6  is  a  constant,  a  design  which  maximizes  e  .  [T-1L(I  8A)L'l 

mm  r 

is  called  a  A^optimal  design.     However,  there  are  some  situations  in 

which    emin[T"1L(Ir«A)L']    is  equal  to  zero  for  any  choice  of  design. 

For  instance  this  happens  when  the  matrix  [Xq:Zq]  is  not  of  full  column 

rank   since  when   [XQ:   ZQ]   is  less  than  full   column   rank    A  will  be 

singular,  and  hence  emin(A)  =  0.     This  implies  that    e  .  (I  8A)  =  0. 

iimii  r  mm  r 

Thus  Ip8A  is  positive  semidefinite,   and  T'1/2L( I  BA)L* (T"1^2) 1  is  also 

positive    semidefinite,    which    implies    that     e  .  [T_1L(I  8A)L'l  =  0. 

mm  r 

Thus,  it  is  clear  that  Aj-optimal  designs  can  be  obtained  only  under 
some  special  conditions.  This  leads  us  to  propose  a  second  design 
criterion  which  can  be  applied  in  more  general  situations. 

3.4.2    Ag -optimal ity 

Our  second  design  criterion  is  to  maximize  A2,  the  average  of  A 
(instead  of  the  infimum  of  A)  over  the  contour  t  =  6,  i.e.,  we  propose 
to  select  a  design  which  maximizes 
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A2  =    J    I'L(IraA)L'IdG/  /  dG  (3.18) 

h 

where  dG  is  the  differential  of  the  area  on  the  surface  of  the  ellipsoid 

^0  =  l'TI  =  5}*  Jones  and  Mitcne11  (1978)  showed  that  the 
following  equality  holds, 

A2  =  q_16A2    where    ^  =  tr[T_1L( Ip8A)L ' ]      .  (3.19) 

A  design  which  maximizes  A2  (respectively  A£)  is  called  a  A2-optimal 
(respectively  A2-optimal)  design.  Since  q  and  5  are  constants  it  is 
clear  that  A2-optimal  designs  and    A2-optimal    designs  are  equivalent. 

If  the  number  of  design  points  N  is  fixed  beforehand  a  A2-optimal 
design  could  be  obtained  by  maximizing  A£  with  respect  to  the  Nk  design 
settings  (coordinates  of  the  N  design  points).  However,  this  may  lead 
to  computational  difficulties  especially  for  large  values  of  N  or  k.  A 
sequential  procedure  by  which  design  points  can  be  chosen  one  at  a  time 
is  therefore  quite  desirable.  Silvey  (1980)  came  up  with  a  sequential 
procedure  to  obtain  a  ^-optimal  design,  if  <j>  is  a  real  valued  function 
of  the  moment  matrix  associated  with  a  linear  single  response  model, 
provided  $  satisfied  certain  properties.  In  our  case  it  can  be  shown 
that  A£  is  a  function  of  the  moment  matrix  associated  with  a  linear 
single  response  model,  and  also  that  AJ,  satisifies  the  required 
properties.  As  shown  in  the  next  section,  Silvey's  procedure  can  thus 
be  employed  to  obtain  a    A'-optimal    design  sequentially. 
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3.5    Sequential  Generation  of  A^-  Designs 
Let  us  consider  the  'artificial'  single  response  model 


£   =  L"x0    Z0]  1  +  £  (3.20) 


where  ya  is  a  Nxl  vector  of  observed  response  values,  XQ  and  ZQ  are  as 
given  before,  _e  is  the  parameter  vector  associated  with  the  model,  and  e 
is  the  Nxl  random  error  vector  with  E(e)  =  _0  and  Var(_e)  =  1.  Let  x  = 
(xj,   x2,  xk)',   where  x_  is  a  design   point   of  the  experimental 

region,  x,  aJQCy)  and  b_'(xj  be  the  row  vectors  which  represent  the  uth 
rows  of  X0  and  ZQ  respectively,  and  H  be  the  set  of  all  design  measures 
on  x  (see  section  2.3.1  for  the  definition  of  a  design  measure).  If 
M(s)  is  the  moment  matrix  associated  with  (3.20),  then  M(?)  is  of  order 
mxm  where  m  =  p  +  pp  and  using  (2.14)  with  r  =  1,  Z  =  1  and  <{>(xj  = 
hjxj,  where  h_' (x_)  =  [a'(x)  :  b'(x)]  we  have 


M(c)  =    E    XJ,(x^)h«(x  )      .  (3.21) 
u=l 

s 

Here,    1  <  s  <  m' ,    0  <  Xy  <  1   with     I    X  =1,    and  m'  =  (m(m+l)/2)+l . 

u=l 

Equation  (3.21)  can  be  rewritten  as 


m1 

M(C)  »  ^  ^Jh'CxJ  (3.22) 

where    X  =  (X     x  \    )'     belongs   to  U  and     U  =  {X:  0  <  X    <  1 

m1^        1     c  n  -  u 

and     z    X    =  1}.    (Observe  that  X    =  0  for    s  <  u  <  m'). 
u  =  l 
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Alternatively, 


M(C) 


Mn(0  Mxz(0 


MZXU)  Mzz(c) 


(3.23) 


m" 


m 


where  Mxx(c)  =  J    ^(x^a '  (xu) ,  Mxz(c)  =    I    ^(xJb'CxJ,  Mxz(c) 

m'  m' 

z    ^(XyJa'Cx  ),        and  Mzz(c)  =    z    KMxJb'U).  Let 

U=l                      u  u  =  l     u-    u  -  -u 

M  =  (M(?):  c  e  H}.  Then  for  M(s)  e  M  we  can  associate  a  point  at  = 
(X^,  x^,  x^,  ,  x^,  x^,,  •••»  2^')'  ^n  the  closed  and  bounded  subspace 
U  x  xm  .  Let  DN  denote  the  set  of  N  design  points  for  fitting  model 
(3.20).  Then  we  can  form  a  design  measure  associated  with  DN,  by 
attaching  probability  1/N  to  each  design  point,  and  using  (2.13)  with 
r=l,  1=1,  and  FD    =  [XQ:  ZQ],  we  have 


M(CD  )  = 
UN 


MXX(?D  >  MXZ(CDJ 


MZX\>  MZZ^DN) 


(3.24) 


where    M^)  -  X^/N,      M^)  -  X^/N,      M^)  -  ZQ'Xn/N, 

Recall   that  A2  =  trCT"1L(IraA)L')],  where  A=Z^(  IN-Xn(X^X0) ^X^Zg. 
Now  since  A  -  N«zz(c,J  -  Nzx(Cn)MjJ(C„  JH^Cc,,  )},     A«  associated 


N 


N 


N 
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with  an  N-point  design  is  a  function  of  the  moment  matrix  M(cn  ).  Thus 
the  generalization   of   the  function   for   any     M(c)  e  M    can  be 

written  as 

AJCM(C)]  =  tP[T-1L{Iri(Mzz(c)  -  Mzx(c)M-J(c)Mxz(c)}L']  . 

(3.26) 

Thus  A£  is  a  function  of  the  moment  matrix  associated  with  a  linear 
single  response  model . 

We  shall  now  describe  Silvey's  sequential  procedure.    The  relevant 

definitions  are  given  first,  and  the  notation  used  is  consistent  with 

model  (3.20).    For  convenience,  we  shall  write    ?n     for  c„. 

N  N 

Let  <j>  be  a  real  valued  function  bounded  above  on  M  but  not 
necessarily  bounded  below,  i.e.,  *[M(e)]  =  for  some  design  measures 
C  e  H.    Let    M'  =  (M(c):  ♦L"M(c)]  >  --}. 

Defini tion  3.1 

♦  is  said  to  be  concave  on  M,  if  for  0  <  a  <  1  and  Ma,  Mb  e  M 
♦  (M)  >  a<f>(Ma)  +  (1  -  a)<j.(Mb),  where  M  =  aMa  +  (1  -  a)Mb. 

Definition  3.2 

Let  M(c)  e  W .  Then  *  is  said  to  be  differentiate  at  M(c),  if  $ 
is  partially  differentiable  with  respect  to  each  element  of  u,  where  u 
is  the  point  in    U  x  xm     which  is  associated  with  M(c). 
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Definition  3.3 

A  design  measure       is  said  to  be  <f>-optimal  if  and  only  if 

♦CM(C*)]  =  sup  *CM(c)]  . 


Definition  3.4 

Let  M^,  M2  e  M.     Then  the  Frechet  derivative  of  <f>  at  Mj  in  the 
direction  of  M2  is: 

F.(M-,  M2)  =  lim+  (l/e)[*{(l  -  e)M.  +  eM?}  -  $(M.)]  . 

For  the  development  of  the  theory  which  leads  to  Silvey's  procedure 
mentioned  earlier,  it  is  necessary  that  (f>  be  concave  on  M  and 
differentiate  on  M'  (Silvey,  1980,  pp.  17-18).  In  addition,  a 
♦-optimal  measure  should  exist  (Silvey,  1980,  p.  22).  The  basic  idea 
used  in  constructing  the  sequential  procedure  (Silvey,  1980,  p.  29)  is 
to  choose  the  design  point  x^+1  such  that 


VM<5N),  h(^+1)h_'(xN+1)}  -  maxF  {M(?  ),  h(x)h'(x)>  . 

XGX  ^ 

The  stopping  rule  is    max  F  {M(c  ) ,  h(x)h ' (x) }  =  0,    which  is  based  on 

xex  * 

the  following  lemma  (Silvey,  1980.  p.  22):  Let  +  be  concave  on  M  and 
differentiate  on  M' .  Suppose  a  <j>-optimal  measure  exists.  Then  c*  is 
♦-optimal  if  and  only  if 
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max  F ,{MU*)  ,  h(x)h'(x)}  =  0 
xex 


The  main  steps  of  the  procedure  are  given  below. 

3.5.1.    The  Sequential  Procedure  to  Obtain  a  <t>-0ptima1  Design 

1.     Start  with  an  initial  design  DNq  such  that    M(cN  )  e  M' 


'0 


2.     Obtain  the  design  point  _Xj\|0+i  at  which 


max  F  [MUN  ),  h(x)h'(x)]  is  attained. 
xcX    9      1  0  ~ 


3.     Obtain   DNq+1    (hence   CNq+1)   by   augmenting  DNq  with  Xf|  +1- 
Recall  that  is  the  design  measure  obtained  by  assigning 

probability  1/(NQ+1)  to  each  design  point  in  DNq+1  . 


4.     Repeat  this  process  to  find  x^0+2,  x^0+3,   ...  by  replacing 


%'  by         ?V2'  "*  until 


max  F  [MU N),  h(x)h'(x)]  <  5 
xex  v 


for  some  N  and  6,  where  6  is  a  small  positive  number  chosen  a 
priori . 
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The  design  measure  Cfj+j.  defined  iteratively  in  the  above 
sequential  procedure,  gives  equal  weights  1/(N+1)  to  each  of  its  support 
points  and  hence  satisfies  the  recursive  formula 

CN+1  =  (1  "  aN)?N  +  V(W 

where     aN  =  1/(N+1)     and  denotes   the  design  measure  which 

assigns  weight  1  to  the  point  Xfl+i*  This  formula  is  interpreted  as 
follows.  cN+1  assigns  probability  aN  to  the  design  point  Xjj+i  and  a 
cumulative  probability   (1  -  oN)  to  DN  =  {xy  %r  ....  x^}  (i.e.,  the 

probabilities  assigned  to    x . ,  i=l,  2,          N    sum  to  (1  -  aN)).  Also 

N            1  n 

11m  aN  =  0    and      I    a     is  divergent  and  for  such  {a.,}    the  above 

N+-    n                  u=l    u  ^ 

sequential    procedure    converges.       This    is    formally  stated    in  the 

following  theorem. 
Theorem  3.2 

The  procedure  described  in  section  3.5.1  converges,  i.e.,  given 
6  >  0, 

max  F  CMUN),  h(x)h'(x)]  <  5 

for  some  N.    See  Silvey  (1980,  pp.  35-36)  for  proof. 

A  well  known  application  of  the  above  procedure  is  Wynn's  algorithm 
(see  Wynn,  1970)  to  obtain  D-optimal  designs.  In  this  situation  * 
sati  sfies 
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log|M(0|    if    M(?)    is  nonsingular 

♦CM(c)]  = 

-°°  otherwise 

Also  F^[M(0,  h(x)h'(x)]  =  h^xjM'^OhJx)  -  m  (see  Silvey,  1980, 
p.  30). 

In  our  application  <|>  satisfies 

A^CMU)]      if  Mxx^)  is  nonsingular 

♦CM(C)]  = 

-°°  otherwi  se 

Hence  M'  =  (M(c):  Mxx(c)  is  nonsingular}.  An  explicit  expression  for 
when  <j>  =  is  derived  in  Theorem  3.3  and  it  is  shown  in  Lemma  3.1 
that  A£  satisfies  the  properties  required  to  employ  the  sequential 
procedure  described  in  subsection  3.5.1.  The  sequential  procedure  to 
obtain  a  AJ,-optimal  design  is  given  in  subsection  3.5.2. 

Theorem  3.3 

Let  ?  e  H.  Then 

FA[M(?),  h(x)h'(x)]  -trtrkc^K^xJ  -  v(x ,  c) )  (b '  (x )  -  v'(x,C)))}L'] 

.  "  A^CM(c)] 

where    v(x,  ?)  =  M, x(c)M"J(c)a(x) . 
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Recall  that    a_'  (x^)  and  b'  (x^)    are  the  row  vectors  which  represent 
the  uth  rows  of  XQ  and  ZQ  respectively,  and    h'{x)  =  [a'(x_):  b*  (x_)]. 
The  proof  is  given  in  Appendix  G. 

Lemma  3.1 

A£  is  concave  on  M,  di fferentiable  on  M'  and  a  A^-optimal  measure 
exists.    The  proof  is  given  in  Appendix  H. 


3.5.2    The  Sequential  Procedure  to  Obtain  a  A^-Optimal  Design 

1.     Start  with  an  initial  design  DNq  (consisting  of  Ng  points) 
for  which  XgXg  is  nonsingular. 


2.     Obtain  the  design  point  _x_Ng+l  at  which 

max  F  ,  [M(cN  ),  h(x)h*(x)]  is  attained. 
xeX    A2  N0 


The  formulas  necessary  to  calculate  F  ,  [M(s    ),  h(x)h'(x)] 

2        N0  ~ 
and  M(cN  )  are  given  in  Theorem  3.3  and  (3.24)  respectively. 


3.  Obtain  DNq+1  (hence  ?Nq+1)  by  augmenting  DNq  with  Xj|  +1. 
Recall  that  cNq+1  is  the  design  measure  obtained  by  assigning 
probability  1/(NQ+1)  to  each  design  point  in  Dn  +1. 


4.     Repeat  this  process  to  find  x^,  x^+3,   ....  by  replacing 

%'  by  ?No+1'  ?No+2'  "*  unti1 


max  F  [M(?  ),  h(x)h'(x)]  <  6 
xeX     2  N 
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for  some  N  and  6,  where  6  is  a  small  positive  number  chosen  a 
priori . 

3.6    A  Numerical  Example 
We   now   give   an   example  to   illustrate  the  sequential  procedure 
described  in  subsection  3.5.1.    In  this  example  we  have  three  responses 
and   two   controllable   variables,    x±   and   x2,   both    in   the   range  of 
[-1,  1].    The  fitted  models  are  given  below. 

yl  =  P10  +  eilxl  +  312x2  +  el 

y2  =  820  +  321X1  +  322x2  +  e2  (3>29) 

y3  =  830  +  031xl  +  B32*2  +  033X1X2  +  e3  ' 

The  true  models  are  all  second  degree  complete  with  pure  quadratic  and 
cross  product  terms  in  Xj  and  x2.  The  matrix  X  is  partitioned  as 
X  =  [Xj:  X2:  X3]  with  both  X^,  X2  having  the  column  of  ones  and  the 
columns  corresponding  to  xj  and  x2  while  X3  has  the  column  of  ones  and 
the  columns  corresponding  to  Xj,  x2,  and  XjX2.  Similarly  the  matrix  Z 
is  partitioned  as  Z  =  [Z^   Z2:   Z3]  with  both  Zj,  Z2  having  columns 

corresponding  to  XjX2,  xj[  and  xjjj  while  Z3  has  the  columns  corresponding 

2  ? 

to  Xj  and  x2  only.     The  matrix  XQ  whose  columns  form  a  basis  for  the 
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columns  of  X  contains  the  column  of  ones  and  the  columns  corresponding 
to  Xj,  x2  anc*  xlx2*    The  matrix  Zg  contains  the  columns  corresponding  to 

op  _ 

x1x2,  x|  and  x£.    At  each  iteration  FA,    is  evaluated  with  the  use  of 

2 

the   expression    derived    in    Theorem    3.3.      The    necessary  quantities 

required  for  the  evaluation  of    F. ,    are  given  below.    The  matrix  T  is 

2 

calculated  using  (3.16).  Since  q  =  8,  r  =  3  and  p1  =  3,  L  =  [Ig:  Ogxl]C 
(see  (3.15)),  where  C  (see  (3.14))  is  calculated  below.  The  matrices 
MZX(?N)  and  MXX(?N)  at  the  Nth  iteration  were  calculated  using 
(3.24).  The  vectors  a_(x_)  and  b_(xj  in  this  example  are  a_'(x_)  =  (1,  X]L, 
x2,  xxx2)  and  b_'  (x)  =  (x^x2,  x\,  x|).  The  matrices  Hi ,  i  =  1 ,  2,  3  (see 
(3.13))  are 

0  0 

1  0 
0  1 

and  the  matrix  C  is  of  order  9x9,  and  is  equal  to  (see  (3.14)) 


and  = 


100000000 
010000000 
001000000 
000100000 
000010000 
000001000 
000000010 
000000001 
000000100 


The  initial  design  consisted  of  the  design  points  (1,  1),  (1,  -l), 
(-1,  1),  (-1,  -1).  For  this  design  XQXQ  is  nonsingular.  The  design 
points  obtained  sequentially  are  given  in  Table  3.1  and  the  value  of  5 
chosen  is  .005. 
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Table  3.1 

The  Augmented  Points  of  the    A^-Optimal  Design 


max  FA,  {M(cN_1),h(x)h'(x)}  A^MU^)] 


5  16.875  (0,  0)  0 

6  8.100  (0,  0)  2.70 

7  3.750  (0,  0)  3.75 

8  1-378  (0,  0)  4.13 

9  -003  (0,  0)  4.22 
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As     seen     from    the    above     Table,      max  FA,  [M(?8),  h(x)h_' (x)]  is 

xeX  2 

approximately  zero.  Thus  we  conclude  that  an  approximate  A^-optimal 
(hence  a  A2~optimal )  design  for  this  example  is  (1,  1),  (1,  -1), 
(-1,  1),  (-1,  -1)  and  five  replicated  observations  at  the  origin. 

As  described  in  section  3.4  the  main  objective  of  constructing  a 
A2-optimal  design  is  to  increase  the  power  of  the  multivariate  lack  of 
fit  test.    This  is  achieved  in  view  of  the  following  facts. 

1.  The  power  is  an  increasing  function  of  the  eigenvalues  of  the 
noncentral ity  parameters  matrix  ft. 

2.  A2-optimal  designs  are  intended  to  increase  the  trace  of  n, 
and  hence  at  least  one  of  its  eigenvalues. 


CHAPTER  FOUR 
DESIGNS  TO  REDUCE  THE  EFFECT  OF  NONNORMALITY 
ON  SOME  TESTS  ASSOCIATED  WITH  REGRESSION  PARAMETERS 

4.1  Introduction 

The  effect  of  nonnormality  on  multi response  regression  tests  was 
studied  by  Mardia  (1971).  Following  Box  and  Watson  (1962),  he  showed 
that  the  null  distribution  of  an  appropriate  test  statistic  could  be 
approximated  by  a  beta  distribution  with  modified  degrees  of  freedom. 
The  only  assumption  made  was  that  the  distribution  of  the  error  matrix 
e'  =  [ej  :  £2  :  ...  :  £fl]rxN  is  symmetric  in  jEj,  e^,  _e^,  where  r  is 
the  number  of  responses  and  N  is  the  number  of  design  settings.  This  is 
true  if  the  error  vectors  are  i.i.d.  Mardia  also  showed  that  the 
modified  degrees  of  freedom  could  be  made  exact  if  a  certain 
multivariate  measure  of  kurtosis  for  the  regression  (controllable) 
variables  (denoted  by  Cx)  is  zero.  Since  Cx  depends  on  the  design 
points  chosen,  it  could  be  made  zero  by  a  suitable  choice  of  design. 

In  this  chapter,  we  investigate  methods  of  design  construction  to 
minimize  |CX|.  This  is  done  by  adapting  the  design  constructing 
techniques  developed  by  Khuri  and  Myers  (1981). 
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4.2    Robustness  Criterion 
Suppose  yu,  u  =  1,  2,  N  is  an  rxl  vector  which  consists  of  the 

observations  made  on  the  r  responses  at  the  utn  design  setting. 
Consider  the  model 

Y  ■  XB  +  e  (4.D 

where  Y  =  [y^  :  y^'  ...  :  yj|]\  X  =  CXX  :  X2  :  ...  :  Xp],  Each  Xi  is  an 
Nxp1-  full  column-rank  matrix:  its  first  column  is  equal  to  l_  (the  vector 
of  ones)  and  the  remaining  elements  are  known  functions  of  the  settings 
of  the  k  controllable  variables,  such  that  the  uth  row  corresponds  to 
the  uth  design  setting.     Also  B  =  diag  (JJj,  j^,  Bp)  where  _B i , 

1*1,  2,  ....  r  is  the  vector  of  regression  coefficients  associated 
with  the  ith  response  and  e  =  [ej,  ^  ej,]'  is  a  Nxr  matrix  of 

errors  such  that  the  vector  e^,  u  =  1,  2,  ....  N  has  the  mean  vector  0 
and  the  variance-covariance  matrix  E,  where  Z  is  unknown. 

Let  the  columns  of  [J_:  0]  form  a  basis  for  the  columns  of  X.  Then 
D  has  order  Nxp  and  rank  p,  where  rank(X)  =  p  +  1.  Also 

Xi  =  H:  D]E1      ,      1-1,  2,  r     ,  (4.2) 


where  Ei  is  a  unique  matrix  of  rank  pi .  In  view  of  (4.2)  we  can  write 
(4.1)  as 
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Y  =  [1:  DXEj  :  E2  :  ...  :  Er] 


i  .e. , 

Y  =  11:  DDCE^:  EgBg :  ...  :  E^]  +  e  . 

Equi valently 

Y  =1  yj+  DG+  e  (4.3) 

where  Yq  is  an  rxl  vector  of  parameters,  G  -  [jfj:  j^:  where  ^  , 

1  ■  1,  2,  r  is  a  pxl  vector  of  parameters  with  some  of  its  elements 

possibly  equal  to  zero.  Let  A  be  a  known  rxm  full  column-rank  matrix 
(r  >  m)  such  that  GA  =  F,  where  F  is  a  pxm  matrix  of  linear  parametric 
functions.    Then  we  can  write  (4.3)  as 

Ya  =  Kjjg)'  +  DF  +  ea  (4>4) 

where  ya  =  YA,  (yjj)'  =  Yq'A  and  ea  =  eA.  It  can  be  easily  seen  that 
each  column  of  (e3)1  has  a  zero  mean  vector  and  a  vari ance-covari ance 
matrix  A'ZA.  Observe  that  the  model  (4.4)  is  of  the  same  form  as  model 
(2.1)  in  Mardia  (1971).    Moreover,  let  M  =  D(D'D)"iD'  and  Z  =  Ya  -1(Y3)' 


h 
0  0 


0 
0 


+  e 
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where  Y3  is  tne  mxl  vector  whose  elements  are  the  column  means  of 
Ya.  Then  an  appropriate  test  statistic  to  test  the  null  hypothesis 
H0  :  F  =  0  against  Ha:  F  *  0  is  ,  V  =  tr(H(E  +  H)"1)  where  H  =  Z'MZ, 
E  =  Z'(I  -  M)Z  and  'tr'  donotes  trace.  Let  s  =  min(r,  p),  P  =  rp/2  and 
Q  =  [s(N  -  1)  -  rp]/2.  Then  Mardia  showed  that,  even  under  nonnormal 
situations  the  null  distribution  of  V/s  is  an  approximate  beta 
distribution  with  the  parameters  6P  and  6Q,  provided  the  joint 
probability  density  function  of    (ea)'  =  [e*:  e|  ejj]    is  symmetric 

in  •••»         Tne  quantity  6  is  called  the  corrective  factor  and 

is  equal  to  one  if  Cx  =  0  where 

Cx=  {N(N  -  1)(N  +  l)/[p(N  -  p  -  1)(N  -  3)]}  x 

(4.5) 

{d  -  [p(p  +  2)(N  -  1)/N(N  +  1)]}  . 

N  2 

In  (4.5),    d  =    z    d. .    and         is  the  itn  diagonal  element  of  M.  Now 
i  =1 

if 


g(D)  =  d  -  p(p  +  2)(N  -  1)/[N(N  +  1)]    ,  (4>6) 

then  g(D)   =  0  implies   Cx  =  0.     In  this  case  6  =  1  and  hence  the 

approximate    distribution    of    V/s    is    identical    to  the  distribution 

obtained  under  the  normality  assumption.     Now  since  d  depends  on  the 

design  points,  an  appropriate  design  will  minimize  |g(D)|.  Such  a 
design  is  said  to  satisfy  the  robustness  criterion. 
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The  quantities  and  g(D)  above  are  identical  to  the  corresponding 
quantities  in  the  single  response  model  (see  Khuri  and  Myers,  1981). 
Therefore,  multi response  designs  which  satisfy  the  robustness  criterion 
can  be  constructed  by  employing  the  same  technique  developed  for  the 
single  response  case  by  Khuri  and  Myers  (1981). 

4 . 3    The  Construction  of  a  Pes ign  Satisfying  the  Robu stness 
Criterion  for  a  Given  Number  of  Experimental  Runs 

In  this  section  we  describe  the  Khuri  and  Myers  technique  for  the 
single  response  case,  and  as  mentioned  earlier,  the  same  technique  works 
for  the  multi response  case. 

Box  and  Watson  (1962)  have  shown  that  d  in  (4.5)  is  invariant  under 

any  nonsingular  transformation  of  the  form  S  =  DT  of  the  columns  of  D. 

Therefore,  the  columns  of  D  can  be  regarded  as  orthogonal.  Suppose 

D  =  [DJ:  (V      ]',    where  D1  is  the  submatrix  of  D  consisting  of  the 
o 

n  =  N  -  nQ  noncenter  points  and  the  controllable  variables  are  scaled 
so  that  D|D1  =  nl.  Let  C  =  [c.jj]  be  an  orthogonal  matrix  of  order  nxn 
whose  first  row  consists  of  elements  equal  to  (1/n)1/2,  and  the  next  p 
rows  (p  <  n)  are  the  corresponding  columns  of  D1/(n)1/'2.  Then 

n     p+1    ?  ? 

g(D)  =    l    (I    C..T  -  p(p  +  2)(N  -  1)/[N(N  +  1)]      .  (4.7) 
j=l    i=2  J 

(See  Khuri  and  Myers  (1981)  for  details.)  In  addition,  the  elements  of 
C  satisfy 
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n     „  n 

Z    c .  .  =  1    ,     Z   c..  »  0    ,    i  =  2,  3,         p  +  1  , 
j=l    1J  j=l  1J 

and  (4.8) 


z    C|jC|j  =  0     ,     2<i<l<p  +  l 


j-1 


Khun'  and  Myers  used  the  following  approach  to  obtain  a  design 
which  minimizes  |g(D)|,  or  equivalently  g2(D),  subject  to  the 
constraints  given  by  (4.8).  Since  g2(D)  is  continuous,  and  the  domain 
of  {c..:  2  <  i  <  p  +  1,  1  <  j  <  n}  satisfying  (4.8)  is  a  closed  and 
bounded  subset  in  the  Euclidean  space,  the  absolute  minimum  of  g2(D) 
must  exist,  and  be  attained  in  that  domain.  If  a  solution  to  g(D)  =  0 
exists,  then  the  absolute  minimum  of  g2(D)  must  be  zero,  and  vice 
versa.  The  method  of  Lagrange  multipliers  can,  therefore,  be  used  as 
follows  to  obtain  a  solution  for  {c.  • :  2  <  i  <  p  +  1,  1  <  j  <  n}, 
which  leads  to  a  design  satisfying  the  robustness  criterion.  First 
consider  the  Lagrange  function 


2         p+1         n     2  p+1  n 


L  =  gt(D)  +    Z     \.(  Z    ct.  -1)  +    E    u.  I 
i=2     1  j=l    1J  i=2 

(4.9) 


P+1  n 
+     2     n.,    z  c..c,. 
i<l=3    11  j-i    ^  ]J 


where  X. ,  y. ,  i  =  2,  3,  . . . ,  p  +  1,  and  n^ ,  2  <  i  <  1  <  p  +  l  are 
Lagrange  multipliers.    Then  solve  the    (2p  +  p(p  -  l)/2)    equations  in 


60 


(4.8),    and   the    np    equations     -,  =  0,     i    =   2,    3   p   +  1; 

9cij 

j  =  1,  2  n  for  c-j  ,•  and  the  Lagrange  multipliers. 


4.4    Finding  an  Optimal  Number  of  Center  Points  to  Construct 
an  Orthogonal  Central  Composite  Design 
Satisfying  the  Robustness  Criterion 

A  central   composite  design   is  obtained  by  adding  the  following 

Is 

design  points  to  a  2  factorial  design,  or  a  fractional  factorial  design 
(see  Myers,  1971.  p.  127). 


xl 
0 

• 

0 

-a 

a 
0 
0 

0 
0 


x2 
0 

• 

0 
0 
0 

-a 

a 
• 

0 
0 


x3 
0 

0 
0 
0 
0 
0 

* 
♦ 

0 
0 


0 

0 
0 
0 
0 
0 
0 

• 

-a 

a 


Thus  if  there  are  nQ  center  points,  then  a  central  composite  design 
which  contains  a  complete  factorial  design,  has  2k  +  n0  +  2k  design 
points. 
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Suppose  the  matrix  D  in  (4.2)  corresponds  to  the  design  matrix  of 
an  orthogonal  central  composite  design,  which  contains  a  complete 
factorial  design.  Then 

g(D)  =  A[(k/(A  +  2a2))  +  (k(l  -  C1)2/*  +  (k(k  -  1)/2A)]2  +  n()k2(C'  )4/<|> 

+  2k[(a2/(A  +  2a2))  +  ((a2  -  C')2/*)  +  ((k  -  1) (C ' )2/*) ]2 
-Cp(p  +  2)(N  -  1)/N(N  +  1)]  . 

Here,  A  =  2k,  *  =  (AT  -  4Aa4)/(A  +  T) ,  C  =  (A  +  2a2)/(A  +  T) ,  p  = 
(3k  +  k2)/2,  and  a  =  (vA/4)1/4  where  T  =  2k  +  nQ,  N  =  T  +  A,  and 
v  =  ((A  +  T)1/2  -  A1//2)2.  Note  that  this  value  of  a  guarantees  the 
orthogonality  of  the  design  (see  Myers,  1971,  Ch.  7).  Now  g(D)  is  a 
function  of  nQ  (  the  number  of  center  points).  Thus  an  optimal  value  of 
nQ  can  be  obtained  by  simply  plotting  g(D)  for  different  values  of  nQ. 
However,  it  is  possible  that  the  g(D)  value  corresponding  to  the  optimal 
value  of  nQ  may  not  be  close  to  zero. 

4 . 5    A  Numeri cal  Example 
This   is   a   very   simple   example  to   illustrate  the   procedure  in 
Section  4.3.     We  have  three  responses  dependent  on  two  controllable 
variables  xj  and  x2.    The  regression  models  are 
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h  -  eiO  +  gllxl  +  P12x2  +  el 

y2  =  320  +  321  xl  +  G2  (4.10) 
y3  =  S30  +  e31  xl  +  e32x2  +  e3  . 

Then  each  of  and  X3  contains  the  column  of  ones  and  the  columns 
corresponding   to     x^,  x^,  contains   the   column   of  ones   and  the 

column  corresponding  to  x..  Therefore  D  in  (4.2)  contains  the  columns 
corresponding  to    x.,  x9    and  p  =  2.    Also  in  (4.2) 


In  (4.3), 


El  =  E3  =  l2     and     E2  = 


1  0 
0  1 
0  0 


G  = 


Pll         821  331 


012        0  632 


We  are  interested  in  testing  the  null  hypothesis 


V    F  =  ?2xl    a9ainst    Ha:    F  *  °2xl 
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where  F  =  (3n  -  g  -  &32) ' .     Then  the  matrix  A  which  satisifes 


F  =  GA  is  A  =  (1,  0,  -1)'.    For  n  =  4,  the  values  of        i  =  2,  3, 


. . , 


p+1;  j  =  1,  2,         n  and  the  Lagrange  multipliers  satisfying    ^     =  0 
and  (4.8)  were  obtained  by  using  the  subroutine  ZSPOW  in  FORTRAN. J  This 
subroutine  solves  a  system  of  nonlinear  equations.    The  second  and  third 
rows  of  C  are  given  by  rows  one  and  two  of  the  matrix 


.500         .500       -.499  -.500 


.499       -.499       .500  -.500 


Therefore 


D'  =  (4)1/2C, 


1  1-1-1 
1-1  1-1 


or    D  = 


1 
1 

-1 
■1 


This  is  a  2L  factorial  design. 


CHAPTER  FIVE 
CONCLUSIONS 


This  dissertation  has  been  primarily  concerned  with  the 
construction  of  designs  for  linear  multiresponse  models.  We  considered 
three  different  topics  in  the  area,  which  are  important  for  better 
design  of  multiresponse  experiments. 

In  Chapter  Two,  D-optimal  designs  for  multiresponse  experiments 
were  discussed,  and  a  sequential  procedure  was  developed  to  obtain 
multiresponse  D-optimal  designs  when  Zt  the  variance-covariance  matrix 
of  the  error  vector  is  not  known.  However,  this  procedure  cannot  be 
used  to  generate  multiresponse  D-optimal  designs  for  a  fixed  number  of 
design  points.  It  would  be  desirable  to  find  a  procedure  that  could  be 
used  to  generate  multiresponse  D-optimal  designs  in  such  situations. 

In  Chapter  Three,  two  design  criteria  were  developed  to  increase 
the  power  of  the  multiresponse  lack  of  fit  test.  They  are  the 
multivariate  extensions  of  the  Aj  and  A2  criteria  proposed  by  Jones  and 
Mitchell  (1978)  for  the  single  response  case.  A  sequential  procedure  to 
obtain  A2-optimal  designs  was  also  presented.  The  Aj  and  A2  criteria 
increase  the  trace  of  q  [the  noncentral ity  parameter  matrix  associated 
with  the  distribution  of  Gj  (see  (3.8))],  and  hence  the  power  of  the 
multiresponse  lack  of  fit  test  which  is  an  increasing  function  of  the 
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eigenvalues  of  fl.  The  question  of  whether  it  is  possible  to  develop 
criteria  to  increase  the  smallest  eigenvalue  of  a  should  be  investigated 
in  the  future. 

In  Chapter  Four  we  considered  the  topic  of  constructing  designs  to 
reduce  the  effect  of  nonnormal ity  of  the  error  distribution  on  tests  of 
hypotheses  associated  with  linear  combinations  of  the  parameter 
vectors.  A  related  topic  for  future  research  is  the  construction  of 
designs  to  reduce  the  effect  of  nonnormal ity  of  the  error  distribution 
on  the  multi response  test  of  lack  of  fit. 


APPENDIX  A 
DEFINITIONS  RELATED  TO  CONVERGENCE  OF 
RANDOM  VECTORS  AND  RANDOM  MATRICES 


Let     (Q,  I,  P)    be  a  probability  space,  and         {x^,  N  >  1}  be 
random  vectors  defined  on  £2  into  R^. 

Definition  A.l 

Xjy  is  said  to  converge  in  probability  to  _x ,  or  equivalents,  is 
a  consistent  estimator  of  x_if  for  any  e  >  0 

lim  P{w:  uiefl  and  I  Ixjju)  -  x(u)||  >  e}  =  0 
N+-  no- 
where   ||  ||    denotes  the  Euclidean  norm  in  Rk. 

Definition  A, 2 

Xfl  is  said  to  converge  with  probability  1  (u.p.l)  to  x_if 

P{u>:  uefl    and    lim    x»  (m)  =  x(w)}  =  1. 
N-m. 

A  matrix  S  of  finite  order  is  said  to  be  a  random  matrix  if  each  of 
its  elements  is  a  random  variable  on  n. 
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Definition  A. 3 

Let  S  and    {SN,  N  >  1}    be  matrices  of  the  same  order.    Then  SN  is 

said  to  be  a  consistent  estimator  of  S  if  and  only  if    s^.     is  a 

■  J 

consistent  estimator  of  s.- •  for  all  1,  j  where  Sw  =  [s^1.]  and 
S  =  &,,]. 

Suppose  S  and  {SN>  N  >  1}  are  symmetric  matrices,  _b(S)  is  the 
vector  which  consists  of  elements  above  or  on  the  diagonal  of  S,  and  if 
eN  is  the  Euclidean  distance  between  b_(SN)  and  b_(S),  then  the 
consistency  of  SN  is  equivalent  to  convergence  of  eN  to  0  in 
probabil ity. 


APPENDIX  B 
RESULTS  IN  MATRIX  ALGEBRA 


Theorem  B.l 

Let  T  be  a  non-singular  matrix  of  order  nxn.  Then 

1.  The  determinant  of  T  denoted    |T|    is  a  polynomial    in  its 
elements. 

2.  Each  element  of  T-1  is  a  rational  function  of  elements  of  T. 

Proof. 

(1)     Using  the  definition  given  in  Aitken  (1951,  p.  30)  we  have 

|T|  =    Z±  V2S        fcnv  (B.l) 

with  the  summation  of  n!  terms  being  extended  over  all  permutations 
(a,  S,  v)  of  column  subscripts  of  the  elements  ti  •  of  T  and  the 

sign  +  or  -  being  prefixed  to  any  term  according  as  the  permutation  is 
even  or  odd,  respectively.    It  is  clear  that  every  term  in  the  summation 
of  (B.l)  is  a  polynomial  in  the  set  of  elements    {t^.,  1  <  i,  j  <n}, 
and  since  a  finite  sum  of  polynomials  is  itself  a  polynomial,  we  have 
the  required  result. 
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(2)     We  have  from  Lancaster  (1969,  p.  36) 


r1  =  (Adj  T)/|T| 

where  (AdjT)^-  =  and         is  the  determinant  of  the  matrix 

obtained  by  deleting  the  jth  row  and  the  ith  column.  Applying  (1)  for 
this  matrix  we  have  that  (AdjT)^  is  a  polynomial  in  elements  of  T  and 
since  |T|  is  a  polynomial  in  elements  of  T  we  obtain  the  required 
result. 


Theorem  B.2 

Let  T  J,  1  <  i  <  j  <  n  be  a  symmetric  matrix  of  order  nxn  such 
that  its  ijth  and  jith  elements  are  equal  to  1  and  all  other  elements 
are  equal  to  zero.  Then  the  non-zero  eigenvalues  of  are  1  and  -1 
for  all    1  <  i  <  j  <  n. 


Proof. 


o  o 


0  0 


0  . 


0  . 


1  j 

0   0 

0    1 

•  • 

1   0 

•  • 

0   0 
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Since  n-2  columns  of  T  are  identically  equal  to  zero,  rank  (T1*  J)  =  2. 
Therefore  T  has  only  two  non  zero  eigenvalues.  Let  g}i  be  a  nxl  vector 
such  that  its  itn  and  jtn  elements  are  equal  to  1  and  all  the  other 
elements  are  equal  to  zero.  Then  £1J  is  an  eigenvector  associated  with 
1.    For  let 


Then 


Therefore 


quj  \lx  'ik  PkJ-  1  <u  <"  where       =  (qJJ    ...    q1J).  , 


qjj  =  tij         =  1 


%J  =  0  for  u  *  i,  j     and  1  <  u  <  n 


Thus    q_ij  =  p1^'. 
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Therefore    T1Jp_1J  =  p_lj  . 

Also  it  is  clear  that  trCPJ)  =  0,  since  all  the  diagonal  elements  of  T 
are  zero.  Thus  if  \  is  the  other  nonzero  eigenvalue  then  1  +  A  =  0, 
which  implies  that  X  =  -1 . 

Theorem  B.3 

Let  S  be  a  pxr  matrix  and  T  be  an  rxr  symmetric  matrix.  Then 
emin(T)  tr(SS'>  *  tr(STS')  <  emax(T)tr(SS' ) 

where  emin  and  emax  denote  the  smallest  and  the  largest  eigenvalues  of  a 

matrix . 

Proof. 


tr(STS-)=trCS(T-emax(T)Ir+emax(T)Ir)S-] 


"  trC-Sta^T)!,.  -  T)S'  ♦  .  <T)SS'] 


<  e  ,  (T)tr(SS') 
max    '    v  ' 


(since  emax(T)I    -  T  is  positive  semidefinite) 
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Simi  larly 


tr(STS-)  -tr[S(T-emjn(T)Ir  tem.n(T)Ir)S'] 


-tr[S(T-e1<1n(T)Ip)S']te1,|n(T)tr(SS-) 


»V,n(T)tr(SS')  . 


Therefore 


emin(T)tr(SS'>  «  tr(STS')  <  emax(T)tr(SS'>  ' 


Corollary  B.l 

Let  A  and  8  be  two  square  matrices  of  the  same  order.  If  A  is 
symmetric  and  B  is  positive  semi  definite,  then 


tr(AB)  <  emax(A)tr(B)  . 


Proof. 

Since  B  is  positive  semidefinite  B  =  S'S  for  some  matrix  S.  Then 
tr(AB)  =  tr(AS'S)  =  tr(SAS')  <  emax  (A)tr(SS')  =  emax(A)tr(B)  using  the 
above  theorem. 


APPENDIX  C  A 
CONSISTENCY  OF  E 


By  definition 


-  (I,  -  WO)  -  Kjl,)     "here  R.  =  F^F'F,)"1^ 

i  ■  1 ,  2 ,  . . . ,  r 


a(£iF;+£i')(IN.Ri)'(IN.R.)(F.Bj+£j) 


-ii(iN-  V(in  -  RjOi.  (C.i) 


N       N  1j 

and  T«  -  [tJJ]  • 
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Therefore 


N 


=  aijtr(T^) 


Now 


tr(T1J)  =  tr[I 


N"Ri  "Rj  +RiRjl 


=  N  -  p.  -  Pj  +  tr(R.R.)  . 

Therefore 

ECN^./CN  -  p.-  Pj  +  trd^Rj)]  =  a.. 


It  is  clear  from  (C.l)  that  7.  is  a  bilinear  form.  Therefore  using 
(52)  in  Searle  (1971,  p.  65)  which  is: 
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Var(x|A12x2)  =  tr(A12C21)2  +  trfA^A^)  +  ^A^C^A^ 

*  u2A21CnA12y2  +  2yJA12C21A12l2  , 

with 

*i  =£i    .   12  -£j    ,    A12  =  (I,  -  R.)(IN  -  R.)  , 
Cll  =  a11Xr    C22  =  ajj!N    •    C12  =  "ij1*  '   ll  =  Jig  =  0 

we  have 

Var[£!(IN  -  R^CIp,  -  R^e.]  -  o*.  trC(IN  -  R,)^  -  Rj)]2 

♦•llV^H-N^II-  V(IH-RJ,(1M-R1>3 

■  4jtr[(IN  -  Ri)dN  -  Rj)]2  ♦  -ii^-jtrCd,  -  Ri)(IN  -  R.)] 
We  also  use  Corollary  B.l  in  Appendix  B  to  obtain 
tr[(IN-V(IN-R.n<eraax(IN-R.>tr(IN-R.)=N.p.  , 


(since  IN  -  R.  is  idempotent,  e       (L.  -  R. 
Mi  max  v  N  i 

Moreover 

tr[(IN  -  R.)(IN  -  R.)]2  =  tr[IN  -  Ri  -  R.  +  R.R.R.R.] 


=  N  -  p.  -  p.  +  tr(R. R.R. R.) 


< N  "pi  -pj  ^(V^V 


=  N  -  p.  -  Pj  +  tr(RjR.Rj) 
Also,  since  (IN  -  R^ )  is  positive  semi  definite, 

RjRj  "  RjRiRj  1s  P°stive  semi  definite. 
Therefore 

tr(R  R  R  )  <  tr(R.R.) 
J  1  J  J  J 


77 

-MR.) 
2 

<  p.      ,     since  p.  >  1 
J  J 

Therefore 

trC(IN  -  R^V  Rj)]2  <  N  -  P.  -  p.  +p2  , 

and 

Var[e'(IN  -  R.)(IN  .  R.)£.]  <  a2.(N  -  p.  -  p.  +  p2)  +  ^..(N  -  p.)  . 
Hence 

Vardifj)  <  o2.(N  -  p.  -  p.  +P2)  +a..a..(N  -p.)  . 
Let  u)sfl  and  e  >  0.    Then  from  the  Chebyshev  inequality 
P{u>:  |Na^(u)/((N  -  p.  -  p    +  trR.R.)  -  a.  .)|  >  e} 

<  Ca2,(N  -  p.  -  p.  +  p2)  +  a..a..(N  .  p.)]/(N  _  p.  .  p_  +  tr(RiRj))2e2  . 
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The  term  on  the  r.h.s.  +  0  as  N  +  <*>.  Therefore  No.  ./(N  -  d.  +  d  + 
trfR^j))  converges  in  probability  to  a^,  1  <  i,  j  <  p.  Also 
(N  -  pi  -  Pj  +  tr(RiRj))/N  converges  to  1,  since  by  Corollary  B.l 
trCR^j)  <  tr(Rj)  =  Pj*  Thus  a.,  converges  in  probability  to  a.j-, 
1  <  1,  J  <  r.    Equivalently    Zw    is  a  consistent  estimator  of  E. 


APPENDIX  D 
PROOF  OF  THEOREM  2.3 


Theorem  2.3 


|M(C,  S)|  =    5    (cr11  )Pi  |M(c,  A"1)],  ;  e  H  and  if1  =  [aij] 
i=l 


Proof.    From  (2.14)  we  have 


M(c,  Z)  =    I    MfeJ  S"V(x,)     where     0  <  X.  <  1  , 
3*1  J  J 


1  <  s  <  p1,  p1  =  p(p  +  l)/2  +1    ,    and    Z    X.  =  1 

Define 


1/2, 


*1  ■  C^^fel).  *2/JW>- .^(5.)]'  p  „.  1  -  1,  2.-.,  r  , 


1/2, 


and  let 


X  = 


X.  0 
X, 


H 
.j 


srxp 
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Then    X'U^BDDX  = 


XJ       0    ...  0 


0  X^  ...  0 
0        0    ...  X' 


0    ...  0 


(2_i8Is) 


0        x2  ...  0 


0        0    ...  X 


orrX'X 
r  r 


j=l    J^   J    1  ~J 


lr 


*     •  • 


s 

Z  X. 

j-1  J 


11 


lr 


^^(XjOfJCXj)  alrfl(x  )f;(x ..) 
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E  X  *(x .)!  $ ' (x.)  ,  where  <fr(x.)  = 
j=l    J     J  -J  -J 


=  M(c,  I),      by  formula  (2.14). 


Therefore 


|M(c,  Z)|  -  |xl (z-1ai)xl 


11 


Let    S  = 


rr 


Then 


A  =  S 


-1/2  r-lc-l/2 


and 


Z-18t    =  S1/2AS1/28I 
s  s 
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=  (s1/2ais)(ABis)(s1/2ais) 


|X'(Z-1SIJX|  =  |X'(S1/28Is)(ABIs)(S1/2BIc)X| 


=  |CX'(ABIc)XC|  where 


C  = 


r(oll)1/ZL 


(a22)1/2! 


Hence 


|  x- (2:-1«is)x  |  =  |C2||X'(ABIS)X| 


r       . .  p. 
■  H    (a11)  1 |X« (A8I  )X| 
i  =1  5 


n   (a1i)Pi|M(?,  A"1; 


APPENDIX  E  A 
CONSISTENCY  OF  A 


Proof 

By  definition 

A  =  (diag  Z"1)'172  E-^diag  z'l)'l/Z 
where    Z  »  [a..]     ,     1  <  1,  j  <  r. 


Suppose    S"    =  [a,J]  then 


Since  a11  >  0,  1  <  i  <  r,  is  well  defined  for  all  1  <  i,  j  <  p. 
Using  Theorem  8.1(2)  we  have  that  a1^  is  a  rational  function  of  an, 
°12'  air  and  nence  eacn  aij-     1  <  i,  j  <  r    is  a  well  defined 

rational  function  of  cru,  o12,  arr.    Thus,  each        is  a  continuous 

function  of  au,  a12,  apr.    Since  is  a  consistent  estimator 
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of  E  and  consistency  is  preserved  under  continuous  transformations,  it 
follows  that    a.,     is  a  consistent  estimator  of  a..--,     1  <  i  <  j  <  r 
where    AN  =  [a..].     Equivalents    AM    is  a  consistent  estimator  of  A. 


APPENDIX  F 
PROOF  OF  THEOREM  2.4 


We  begin  by  introducing  some  notation  and  definitions  related  to 
Theorem  2.4  and  then  establish  the  lemmas  that  are  required  for  its 
proof. 

F.l  Notation 

Let  A(Z)  be  the  set  of  all  matrices  of  the  form  M(c,  £)  where  ?eH, 
the  set  of  design  measures  defined  on  X'    Then  from  (2.14)  we  have 


(F.l) 


s 


where  0  <  s  <  (p(p  +  l)/2)  +  1 


0  <  Xu  <  1,    and    £    X   »  1 


This  can  be  rewritten  as 


(F.2) 


where    p'  =  (p(p  +  l)/2)  +1     and     0  <  X    <  1     with     Z    X    =  1 

u  ,  u 

u=l 

(Observe  that  Xu  =  0  for    s  <  u  <  p'). 
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Let  A  be  as  in  Section  2.6  and  let  A'  be  the  subset  of  A  which 
consists  of  all  positive  definite  matrices  in  A.  Then  A(A)  is  defined 
to  be  the  set  of  all  matrices  of  the  form  M(c,  A-1),  AeA'  given  by 

-1  P' 

M(c,  A    )  =    Z     X  *(x  )A*'(x  )  (F.3) 
u=1     u    -u  -u 

P' 

where    0  <  X    <  1,  and    z    X  =1. 

u-1  U 

Observe  that  (F.3)  was  obtained  on  replacing  Z"1  by  A  in  (F.2),  and 
that  A(A)  is  a  convex  set. 

Let  {Xjj,  N  >  1}  be  a  sequence  of  design  points  in  x  and  U^,  N  >  1} 
be  a  sequence  of  discrete  designs  such  that  the  spectrum  of         is  the 
union  of  the  spectrum  of  eN  and  x^.     If  N'  >  1  and  CN«  is  a  design 
such  that  M(cN,,  A'1)  e  A(A)  is  positive  definite,  then  for    0  <       <  1 
and  N  >  N',  we  define 

\<W  A_1)  =  U  -  A_1)  +  V^W^W'  N  >  N' 

(F.4) 

VaM    <*'  CN'  A"!)  *  ♦'teX1    (5N,  A-1)*(x),    x  e  x.  N  >  N'  . 
N-l  ™  —  — 
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Lemma  F.l 

(1)  {Ma  (5  N+1>  A    ),  N  >  N'  -  1}    is  a  sequence  of  positive  definite 

N 

matrices  in  A(A) . 

(2)  The  elements  of  Ma  A"1),  N  >  N1  -  1  are  polynomials  in 
elements  of  _b(A),  the  r' -dimensional  vector  associated  with  A 
where  r'  =  r(r-l)/2. 

(3)  |Mo  (CN+1.  A"1) |,  N  >  N'  -  1  is  a  polynomial  in  the  elements  of 
b_(A). 

(4)  The  elements  of  M~*UN+1.  A"^'  N  >  N'  -  1  is  a  rational 
function  in  the  elements  of  b(A). 


Proof 

(1)     We  shall  prove  the  result  by  the  method  of  mathematical  induction. 
For  N  =  N'  we  have 


N-l  N 

Observe  that  by  definition    M(cN',  A"1)    is  positive  definite  and 

belongs    to   A(A) .      Suppose   that     M       (5     A"1)     is  positive 

Vl  n 

definite  and  belongs  to  A(A).    We  also  have 


MaN^N+l'A_1)  "  (1  -  ^\_^  A"^  +  V^IM  )A* '  ^N+l )  • 

and  since  M  (cN,  A"1),  ♦(xN+1)A*'(xN+1)  e  A(A),  the  convexity 
of    A( A)    implies  that    M    (eN+1>  A"1)  e  A(A).      Next  we  need  to 
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show  that  Ma  UN+1>  A-*)  is  positive  definite.  This  is  done  by 
means  of  the  following  result  in  Graybill  (1969,  p.  330): 

If  P  is  positive  definite  and  Q  is  positive  semidef inite  then 
|P  -  Q|  <  |P|. 

In  our  case,  P  =  M    (cN+1,  A"1)  and  Q  =  aN*(_xN+1)A*'  (xj,+1) 
and  we  have 


N  N-l 

Now  since    M       (cNt  A"1)    is  assumed  to  be  positive  definite  and 
N-l  n 

«N  <  1  we  have     IM^U^,  A_i)|  >  0.    Since    M    (?N+1,  A'1)  is 

positive  semidefinite  we  conclude  that    M    (cu      A"1)    must  be 

aN 


positive  definite. 


Let  N  >  N'  -  1.    We  have 


P' 


where    0  <  XHt+1  <  1  with    Z    XN+1  =  1. 


u  i  u 

u=l 


Also 
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5.  I2<>T> 


.  ,  N+l. 


1      3  i  0    .  .  .  3  , 

12  lr 


a21    1  ••*a2r 


lrl 


l}(x^+1)     0'  ...  0' 


0' 


arlrr^u    ,fl(*u  ) 


alr^l(^j    'V^  } 


(3) 
(4) 


It  is  clear  from  the  above  matrix  equality  that  each  element  of 
N+l  N+l 

^\    )A+'(Xy    )    is  a  polynomial  in  a12,  a13,  ar_1)P  which 

are  the  elements  of  b_(A) .    Now  from  (F.5)  it  is  clear  that  each 
element  of  M    cN+1,  A-1)  is  a  finite  linear  combination  of  the 
corresponding  elements  of     *(x^+1)A<>'  (x^+1) ,    u  =  1,  2,  p' 
and  hence  is  a  polynomial  in  elements  of  _b(A). 
Follows  from  (2)  above  and  Theorem  B .1(1) . 
Follows  from  (2)  above  and  Theorem  B.l(2). 
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Lemma  F.2 

Let  ceH  amd  AeA1 .    Then  each  element  of  M(c,  A"1)  is  a  polynomial 

function  in  xj,         ...  x,,.,  \v  X2,  xp,,  a12,  a13,    ....  ap_l  r 

defined    on      XP  x  U  x  [-1,  l]r  ,     where    p«     =    r(r-l)/2,    and    U  = 

p' 

IX:  X  =         x2,  ...  x  ,)  and     0  <       <  1    with     z    X    =  1}.  • 

u=l  u 

Proof 

By  (F.3)  we  have 

-1  p' 
M(c,  A    )  =    E  ^♦(xu)A*'(x) 
u=1    u   -u  -u 

P' 

where    0  <  X    <  1  with    I    X    =  1. 

u=l  u 

As  shown  in  the  proof  of  Lemma  F.l(2)  we  have 
♦  (xJAVCxJ 


a  ,f  (x  )fi(x  ) 
r  1-r v-u  -lv-u ' 


f  (x  )f  (x  ) 
-r  -u  '-r  v-tr 
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Therefore 


M(C,  A"1) 


r 


U -i  U  =1 


p' 


p' 


U-l 


J 


Now  f.j(x_,),  i  -  1,  2,  r  is  a  polynomial  function  in  x_  and  hence  it 

is  clear  from  the  above  equality  that  each  element  of  M(c,  A"1)  is  a 
polynomial   function  in         Xg,  x^,  X1§  X2,  xp,,  a12,  a13, 

ar_i>r  defined  on    xP   x  U  x  [-1,  l]r'. 

Lemma  F.3 

Let  ceH  be  such  that  M(c,  A51)  is  positive  definite  for  some 
AQe  A'.     Then  M(?,  A"1)  is  positive  definite  for  all  A  e  A' . 


Proof. 


M(C,  A"1)  -  ^  Xu*(xu)A0^(xu) 


P' 


where  0  <  X    <  1    with     Z    X    =  1 . 
U  u=l  u 
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We  can  also  prove  that 


F'  (A  8  Ip,)F 


where  F  is  a  rp'xp  matrix  such  that 


F  = 


F.       0    ...  0 


*  F 


in  which  is  a  p'xpi  matrix  of  full   column  rank  whose  u 


th 


row 


contains  the  elements  of    ^^(x^),    u  -  1,  2  p'.  Thus 


H(C,  A'1)  =  F'(ABI)  F 


-  F'(Aj/29Ipl)(Aj/2XIpl)F 


1/2 

(AQ'    exists  since  AQ  is 
positive  definite) 
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Also  if  AeA'  then  A1/2  exists  (since  A  is  positive  definite)  and 
rank[(A1/2BIp.)F]  =  rank[(A1/2BIp, ) ( A"1/2XIp . ^( A^7 2XIp , )F] 

=  rank[(Aj/2BIp,)F] 

=  rank[F'(A0BIp,)F] 

=  rank[M(c,  A"1)] 

-  p    (since,  M(c,  A"1)  is  positive  definite). 

It  follows  that 

rank[M(c,  A'1)]  =  rank[F' (All  , )F] 

=  rank[(A1/2gIpI)r] 
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Now  since  M(?,  A"1)  is  of  the  form     F'  (A1/2BI)(A1/28I)F,     it  is 
positive  semi  definite  and  therefore  positive  definite. 

Lemma  F.4 

Let    AQ  e  A'    and  cN'  be  a  discrete  design  such  that    M(?N,,  A"1) 
is  positive  definite  and  let    {AN>  N  >  N ' }    be  any  sequence  in  A'  such 
that  _b(AN)  converges  to  b_(AQ) .     Then  for  any  aN  such  that  0<      <  1, 
there  exists  a  positive  integer  Nj  such  that 

^\^n+v  A',1)  |  >  loglM^C^,  A-^l  + 

VtrCVVl0w  V  an^] 

-p}  +  rN  -  2r'pe-i1n  (AQ) (e°N  +  e^) ,  N  >  H 


where  rN  =  0(oN),  e°  denotes  the  r' -dimensional  Euclidean  distance 
between  _b(AN)  and  _b(A0),  and  emin  denotes  the  smallest  eigenvalue  of  the 
matrix  inside  parantheses. 

Proof 

We  shall  first  prove  that   M  A'1)    is  positive  definite  for 

N-1  n 

any     b(A)  e  B(A     6)     and     N  >  N      where    B(AQ,    5)    is    a  spherical 
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neighborhood  of  radius  <5  centered  at  _b_(A0),  and  N-,  is  a  positive 
integer.  Now  AQ  is  positive  definite  implies  emin(A0)  >  0,  and  since 
emin(A)  is  continuous  in  the  elements  of  A  there  exists  5  >  0  such  that 

emin(A)  >  emin(V/2  >  0  (hence  A  is  positive  definite)  for  all  A,  such 
that  b_(A)  e  B(AQ,  6).  Therefore,  A  e  A',  whenever  _b( A)  e  B(AQ,  6). 
Also  since  b_(AN)  converges  to  b_(AQ)  there  exists  N1  >  N'  such  that 
b_(AN)  eB(Aq,  6)  for  all     N  >  N^.      Now  since  M(cN,,  A'1)    is  positive 

definite   by   hypothesis,    it   follows   that     M       Uw,  A'1),  N  >  N1  is 

aN-l 

positive  definite  in  A(AQ)  in  view  of  Lemma  F.l(l).  Moreover,  since 
b_(A)  e  B(AQ,  5)  implies  AeA'  we  have  from  Lemma  F.3  that  M  A"1), 
N  >  N1  is  positive  definite  whenever  _b( A) eB(A0,  6).    This  implies  that 

l09'MaN-l^N'    A~^'    iS   Wel1    defined   wnenever   b_(A)    e   B(A0,    5)  and 
N  >  Hy    Let    N  >  H1  +  1    be  fixed.    Then  Jb(AN),  b/A^)  e  B(Aq,  5)  and 
for  fixed  oN-1,  log|MaN    (cN,  A_1)|   is  a  polynomial   function  of  _b(A) 

defined  on  B(AQ,  5).    Hence  we  denote  log|M_     (?N,  A_1)|  by  h(A).  Now 

N-l 

h(A)  is  differentiable  on  B(AQ,  6)  and  therefore  we  use  the  Mean  Value 
Theorem  in  many  variables  (see  Gillespie,  1951,  p.  60)  to  write 

h(A)  =  h(A  )  +   Z    (b  -  b")[JU(A)]  ,    where  0  <  8  <  1, 

l=l     1     1    dbi  a=A(9) 

b(A)  e  B(AQ,  5)    with    b(A0)  =  (b°,  b°,  b°.)'  (F.6) 

and  b(A(9))  =  (1  -  e)b(A)  +  9b(A  )  . 
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Notice  that  _b(A(e))  e  B(AQ,  6).    Next  we  prove  that 


g  1 

|yb-h(A)|  <  2pe^in(AQ)  whenever  b(A)  e  B(A  ,  6),  and     1  <  1  <  r\  Now 


,5^(A)  -,Bjiog|H       UN,  A"1; 


"l"Vl(!  N-  A      l-Sbfa,.^  IT  (F.7) 
(using  equation  1.134  on  p.  21  in  Fedorov  (1972)). 


Also  from  Lemma  F.l(l)  we  have 

MViUn'  A"1}  =  u-i  XI*fcU)A*'(*U)     With  0  <  X!  <  1  and 


Therefore 


P'  N 

z    xn  =  1 

u-1  u 
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=    I  X 


N  3 


ll-l    U  3bi 


u=l 


u  mj  '    r  v-u 


,  where  A=[a1-j] 


where  B*  is  the  rxr  symmetric  matrix  such  that  its  ijth  and  jith 
elements  are  equal  to  1  and  all  other  elements  are  equal  to  zero,  since 


.-1- 


bA  is  equal  to  some  a^..    Substituting  for   ^~  M  A"1)    in  (F.7) 


we  obtain 


^  h(A)  -  trCM^U,,  A"1)  ^       (Cn,  A"1)}] 


P' 


=  trCM"1    (cN,  A-!){  r  A(xN)bV0<N)}] 


vr'N 


u=l 


U    X-U '      r  XHJ 


tp[M,M'{  I    A(xlJ)B%'(xlJ)}.  where    M"1    (?w,  A-1)  = 
u«l  ~^  aN-i  N 


MiMi 
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=   I    Xutr[M^(xuN)B%'(xuN)M1]  . 

Hence 

p  i 

I  u — 1 

<  j   ^|trCM*«(xJ)B%«  (x^JMjDl     .  (F.9) 

Now 

*  •«ix(B*)trCMi*feJ)*,0i),,i3  •  i  < "  <  p' 

(using  Theorem  B.3) . 

Thus 

-trCM}*^)*'^]  <  trCM^^B^'^M^  <  tr[K{*(xJ)*'  (xJJjMj] 
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(since  emin(B*)  =  -1  and  emax(B*)  =  1  from  Theorem  B.2), 

and  so 

|tr[M^(xU)B%'(x|J)M1]|  <  trCMJ*(^)V(xuN)M1]  .  (F.10) 
Moreover 

emin(A)trCMi*(xJ|)V(xuN)M1]  <  tp[M^(xJJ)A*' (xjJjMj] 

(using  Theorem  B.3). 

It  follows  that 

trCM^V'^^]  <  e-J^AJtrCM^^Af^jM  J 

<  2e;i1n(A0)t^Mi*^N)A*'^N)Ml^ 
(recall  that  emin(A)  >  emin(A0)/2  >  0,  whenever _b(A)  e  B(Aq  a)). 
Now  using  inequality  (F.10)  we  get 


Finally  using  the  above  inequality  together  with  (F.9)  we  obtain 

h*-h(A)| 

N-l  u=l  ^ 

"  2e™VA0>tr^    (CM.  A-1)Ma      (C  A"1)] 
N-l  N-l 
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Consequently 


|g£-h(A)|  <  2e~11n(A0)p  whenever  b(A)  e  B(AQ,  6),  and 

1  <  %  <  r'  . 

In  particular,  |[-i-h(A)]  |<  2e~?(An)p  .since  b(A(9))  e  B(A,  6) 

°l  A=A(8)        mln  0 

Now  from  (F.6)  we  have 


h(A)  =  h(AQ)  +    z    (bt-bjj)  g^h(A(9)) 

where  Cgl-h(A)]  =  J-  h(A(8)) 

dD*        A=A(9)  3b£ 


Thus 


|h(A)  -  h(AQ)|  =  |  E        -  bj)  ^-h(A(6))| 


*  *.  i">»-i>JihK*(*(«))i 

*=1  % 


r  ■ 

<    z  2e~Jn(AQ)p  e°  where  e°  is  the 


Euclidean  distance  between  b_(AQ)  and  _b(A) 
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|h(A)  -  h(AQ)  |  <  2peVjn(A0)r'  whenever  b_( A)  e  B(AQ>  6). 

Recall  that  b_(AN),  bjA^)  e  B(AQ,  6).  Therefore,  replacing  A  by  AN  and 
AN-1  we  9et 

IMV  -  h(A0)|  <  2pefo11n(A0>r'  , 
and 

'tl'AN-l'  '  h'*o'l  *  2peS-lem1n'*o'r'  • 
These  two  inequalities  imply  that 

-2"r'ek"!n<Ao)  *  h<V  "  h(A0>  *  ^folVV  > 
and 

-2pr'eN-le™"lVA0»  *  h<Vl>  "  h'V  <  2P'-'^.1e„-11„(A0) 
from  which  we  obtain 

h(AN)  -h(AN_l}  >  -Zpr^VV^^S-i)  ' 
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Substituting  for  h(AN)  and  h(AN_1),  we  finally  have 

]°*\Jh>  AN^I  >  l0glMaN.lUN'  ^  '  "^'^  V (*"  +  ^ ' 

N  >  N1  +  1    .  (F.ll) 

The  rest  of  the  proof  somewhat  parallels  the  arguments  given  in  Silvey 
(1980,  p.  36): 

The  matrix    Ma  (?N+1»  Aj^1)    was  defined  as 
N 

For  given  N  >  Nj  +  1,  log|M  (?N+1,  A'1) I  is  a  di fferentiable  function 
of  aN  and  therefore  we  can  use  Taylor's  Theorem  about  oN  =  0  to  get 


«H^°9|MaN(CN+1.  A'1)!]  +rN 

V° 

where  rN  =  0(«N) 
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(using  equation  1.1.34  on  p.  21  in  Fedorov  (1972)) 
l0g|MVl(CN'  AN^I  +  Vtr^_^N'  AN1)*^+l)V'^+l)>P> 

l09,Vi(C""        +aNitr{V'^i)M;i/r  aS1)*^+i)3-p1 

^'Vl^'  AN-l)|-2Pr,(eS+eS-l)^n(V  + 

VtrCVVl<>W  CN>  AN^  "  Pi 

+  rN,    (using  inequality  F .11) . 

We  are  now  ready  to  prove  the  main  result. 
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Theorem  2.4 

N  . 

Suppose    that      z  e      converges    in    probability    to   some  random 
u=l 

variable  e,  then  for  a  given  5  >  0,  there  exists  an  integer  N  >  0  such 
that 


sup  tr[ANV(x,  cN,  A"1)]  -  p  <  6 
xeX 

with  probability  1  (see  (2.20)  for  the  definition  of  AN;  Section  2.6 
for  the  definitions  of    ?.„  e.,) . 

IN  N 


Proof 

N 


Since     ey  >  0,    I    e      is    a    nondecreasing    sequence    of  random 


u=l 
N  . 

variables.     Thus     {  Z    e     N  >  l}    is  a  nondecreasing  random  sequence 

u=l 

which  converges  in  probability  to  a  random  variable  e,  and  hence    £  e 

u=1  u 

converges  w.p.l  to  the  same  random  variable  e  (Bil lingsley,   1979,  p. 

N  . 

236),   i.e.,   if    QQ  =  {u:  ijm    E    e  (u)  =  e(u)j>     then  p(a  j  a  u  Lftt 

»  e  flQ.     Then     Tim  eN(«)  =  0.     Equivalently ,     Mm  b(AN(o>))  =  b(A). 
N+-  [\|x» 

Using  Lemma  F.4  with  AN  =  AN(a,),  aN  =  l/(N+l),  CN  =  ?N(u>),  and  Ag  = 
A  (as  defined  in  (2.19)),  we  have 
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log|M(cN+1(«),  A-^a.))! 


>  log|M(cN(»),  A"1N_1(«)) |  +  (l/N+l){tr[AN(w) 


V(2Lm+i(^)  »  CN(«),  A~a(co))D  -  p} 

+  rN  "  2r'P  %rfn<A,<V">  +  ViW)    •    N  >  V-) 
(note  that  when  aN  -  1/(N+1),  ry^.A"1)  and  V^fx,  cN,  A"1)  reduce 

to  M(cN,  A"1)  and  V(x_,  ?N,  A"1)  mentioned  in  Sections  2.3-2.6) 

log|MUN(a>),  A"^(»))|  +  (l/N+DCsup  tr{AN(o,)V(x,L((o),A:1(a,))}  -  p] 

_xex  n 

+  rN  "  2r'P  e~Jn^^eH^  +  Vj^))    ,    N  >  NjM  . 


Recall  that    trCA^Qc^ ,  C|(.  A"1)]  =  sup  tr[ANV(x.  L.  A"}] 

X£X 

(see  Section  2.6). 

The  proof  from  this  point  on,  closely  follows  the  arguments  given  in 
Silvey  (1980,  pp.  35-36).  Suppose 
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sup  tr[AN(o»)V(x,  cN(«»)f  A'^oi))]  -  p  >  5  for  all  N  . 
xex 


Then 


log|M(JN+1(»),  AjJ^o.))! 


>  log|M(JN(»),  A-^H)!  +  5/(N+l)  +  rN-  2r'pe^n(A)  (eN(<o)  ♦  e^to)), 


N  >  N1(o>)  . 

Using  the  above  recursion  formula,  and  by  repeated  substitution  we  have 


log|M(CN+1(*))f  A-1(a))|-log|M(cNi(uj)+1(a)),  ^      („) )  | 


-2r'pe~|  (A)(    z  J  (.,  +   !         *  (•))  + 

u=N1(o))+l       u      u-N1(u)  u 


N  H 

62  1/(U  +  1)  +    z  r      where     p    =  0(l/u+l) 

u=N1(u>)+l  u=N1(o)+l    u  u  ' 


u=N1(u)+l  u=N1(<u) 
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N  N  N 

*    \  t  *  /1/u+1){6  +  (  1    ,         r J  Z  (1/u+l))}  . 

u=N1(a))+l  U=N  (M)+l    u  u=N1(u))+l 


The  limit  of  the  second  term  on  the  right  hand  side  as    N+«»   is  «  (see 

Silvey,  1980,  p.  36).     Therefore,  the  limit  of  the  right  hand  side  as 

N  ♦  »  is  -2r'pe[;.1n(A){2e(a,)  -  2  I  e      i    .    }  +  .  .        impiying  that 

u=l  i v  ^ 

*— 1 

log|M(?N+1(«u),  A"  (u))|     diverges   to  »  for   any     we^.  Equivalently 
*— 1 

1og|M(cN+1.  A~  )|  diverges  to  -  w.p.l.  But 

log|M(c,  A    )| ,  C  e  H,  A  e  A'     is  a  continuous  function  on  the  closed 
and  bounded  subspace    XP"  x  U  x  [-1,  l]r'.      Hence    1og|M(cN+1,  A"1)! 
is    a    bounded    random    function.       This    contradicts    the    fact  that 
1og|M(?N+1,  A~  )|    diverges  to  -  w.p.l.    This  contradiction  establishes 
the  fact  that    sup  tr[ANV(x,  ?N,  A"1]  -  p  <  6    for  some  integer  N  w.p.l. 


APPENDIX  G 
AN  EXPRESSION  FOR  FA,  [M(c),  h(x)h'(x)] 


Only    an   outline   of   the   proof    is    given   and    routine  algebraic 
manipulations  are  omitted. 
By  definition  3.4 


F  ,  (M,  hh1)  «  lim,  (1/e)  [A'(M)  -  Ai(M)] 
2  e-K)  c  L 


where  M  =  (1  -  e)M  +  eh  _h' . 
Recall  that 


AJ(M)  =  trCT1WIrKlla  -  HZXM>XZ»L'] 


Therefore 


F  ,(M,  hh1)  =  11m,  trCT^Kl  BE)L'} 
2  e+0  r 


where  E  =  -  -  Mzz  +  } 
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and 


M  M 


Mzx  Mzz 


Using  the  following  identity  from  Dykstra  (1971)  for  a  nonsingular 
matrix  A: 

(A  +  XnxA)-1  =  A"1  -  A 

l+x'Ax 

with 

A  =  (1  -  e)Mxx     and    x^  =  e1/2a 

we  have 

M"J  =  [(1  -e)M     +  eaa']"1 


=  (1  -  -  Uc/(1  -  e)2}M-Jaa'M-J 


where  c"1  =  1  +  {e/(l  -  e)}a'M^a 


Ill 


i  .e. , 


*xx  =  (1  "  e)"1mxx  "  {e/t(1  "  e)}P 


where  t  =  1  -  e  +  ea'M~xa_    and    P  =  M'^aa'M"*  . 


With  the  use  of  the  above  equality  and  cancellation  of  terms  we 


E  =  -Mzz  +  MzxMx>xz  +  ^  -  MzxMxx^-'  -  ^'Mx>xz 


-{e/(l  -  e)}ba'M^ab'  +  {(1  -  e)/t}MzxPMxz 


+  (e/t)MzxPab'  +  (e/t)ba'PMxz  +{e2/(l  -  e)t}ba'Pab' 


Using  the  above  equality  and  evaluating  the  limit  of  each  term 
tr[T-1L(I  K)L']    as  e  goes  to  0+  we  obtain 


FA,[MU),  h(x)h'(x)]  =  tr[T"1L{Ir8((b(x)  -  v(x,  c)) 


(b'(x)  -  v'(x.  0))>L']  -  Ai[M(c)] 


where  v(x,  C)  =  Mzx(c)M^(c)a(x) . 


APPENDIX  H 
PROOF  OF  LEMMA  3.1 


Lemma  3.1 

A^    is  concave  on  M,  di fferentiable  on  M1 ,  and  a  A^-optimal 
measure  exists. 


Proof 


Recall  that 


M  =  (M(c),  C  e  H},  and  M'  =  {M(c):  Mxx(c)  is  nonsingular}  . 


Suppose  that    »  =  (X^  X2,  Xm,  ,  x^,  x.2,  ....  x^,  )'    is  the  point  in 


U  x  xm     which  is  associated  with  M(c)  e  M.  Then 


N(C)  ■ 


MXX(C) 


MXZ(C) 


_  Mzx(?)  Mzz{5) 

where  Mxx(0-    =    V^Ja'tx^),  M^(c)  -    E    ^a(x^ )b 1  (x^ ) .  MZX(C) 
u=l  u=l 
m'  m' 

^(xja'^),  and    MZZ(C)  -    X  )b ' ) . 

Also 


AJCM(5)]  =  trCrkd^M^tO  -  Mzx(c)MxJ(5)Mxz(c))>L']  . 
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It  is  clear  that  the  elements  of  Mxx(c),  Mxz(0.  Mzx(0,  and  Mzz(?)  are 
polynomial  functions  of  u  and  hence  A£[M(0]  is  a  continuous 
function  on  the  closed  and  bounded  subspace  Ux  x"'  •  Hence  a 
A£-optimal  measure  exists,  since  a  continuous  function  on  a  closed  and 
bounded  subspace  attains  its  optimum  value. 

Next,  the  differentiability  of  A£  on  M'  is  proved  as  follows.  We 
know  that  the  elements  of  Mxx(?),  MXZU),  Mzx(s),  and  Mzz(c)  are 
polynomial  functions  of  w.  Now  from  Theorem  B.l(2)  we  also  know  that  the 
elements  of  Mxx(c)  are  rational  functions  of  the  elements  of  Mxx(0. 
These  facts  imply  that  the  elements  of  the  matrix  T-1L{I  ■A[M(c)])L,t 
where  A[M(c)]  =  Mzz(?)  -  MZX(OM~xU)MxzU)  are  rat1onal  functions  of 
w.  Since  A£[M(c)]  -  trCT^U^iACMCcJDL' ]  and  the  trace  of  a  matrix 
is  the  sum  of  its  diagonal  elements  A^[M(c)]  has  partial  derivatives 
with  respect  to  all  the  elements  of  oj. 

We  shall  now  prove  that  A£  is  concave  on  M,  i.e.,  for  0  <  a  <  1 
and  Ma,  Mb  e  M 

A^(M)  >  (1  -  a)A^(Ma)  +  aA^(Mb) 

where  M  =  (1  -  a)Ma  +  oMb.  We  shall  use  the  following  inequality  from 
Fedorov  (1972,  p.  20): 

If  Aj  has  dimension  nxm  and  Bj  is  a  square  positive  definite  matrix 
of  order  mxm  (j  =  1,  2)  then 
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(1  -  a)A1B[1A^  +  aA^"^ 


>  [(1  -  o)Aj  +  aA2][(l  -  a)B1  +  aBgrtd  -  <*)A^  +  aA^]  , 

0  <  a  <  1      .  (G.l) 


We  have  that 


MZZ  "  MZXMXXMXZ 


=  (1  -  a)M^  +  aM^  -  [(1  -  a)M^  +  aM^] 


[(1  -  a)M^x  +  aM^rtd  -  «)mJz+  M^] 


>  (1  -  a)Mazz  +  aMbzz  -  (1  -  a)Mazx(Maxx)-Vxz  -  oM^M^)  ^ 


(applying  inequality  (G.l),  with    ^  =  Mzx,  A2  =  Mzx, 


Bl  =  MXX'  B2  =  MXX) 


(1  -  a)[M*zz  -  Mazx(Maxx)-Vxz3  +  «CMbzz  -  ^Ax^ul  • 


Hence 
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>  (1  -  «)trCT"1LCIpi(H|2  -  MzX(MXX)_1mxZ)}L,] 

+  a  trCT-kd^  -Mbzx(Mbxx)-1M5z)}L'] 
This  inequality  implies  that 

AJCM]  >  (1  -  a)A^[Ma]  +  oA^[Mb]  . 
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